【編集履歴】

  • 2023.08.16:執筆開始
  • 2023.08.16:「挿入ソート」を追加
  • 2023.08.20:「マージソート」を追加
  • 2023.08.30:「バケットソート」を追加
  • 2023.09.02:「クイックソート」を追加
  • 2023.09.28:「ヒープソート」を追加

12.3 挿入ソートの実装と可視化

 挿入ソート(insertion sort)を実装します。また、ソートの推移のアニメーションを作成して、アルゴリズムを確認します。

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(gganimate)

 この記事では、基本的に パッケージ名::関数名() の記法を使っているので、パッケージを読み込む必要はありません。ただし、作図コードについてはパッケージ名を省略しているので、ggplot2 を読み込む必要があります。
 また、ネイティブパイプ演算子 |> を使っています。magrittr パッケージのパイプ演算子 %>% に置き換えても処理できますが、その場合は magrittr を読み込む必要があります。

 「ソートの可視化」において、gganimate パッケージを利用してアニメーション(gif画像)を作成します。ソートの実装自体には使いません。

ソートの実装

 まずは、挿入ソートのアルゴリズムを関数として実装します。

 挿入ソートを実装します。

# 挿入ソートの実装
insertion_sort <- function(vec) {
  
  # 数列の要素数を取得
  N <- length(vec)
  
  # 要素ごとに処理
  for(i in 2:N) { # (1番目は不要)
    
    # i番目の値を取得
    val <- vec[i]
    
    # i番目から逆順に挿入位置jを探索:(j ≤ i)
    for(j in (i-1):0) { # (i番目は不要・0番目は挿入処理との兼ね合い用)
      
      # j番目とi番目の値を大小比較
      if(j <= 0) {
        
        # 全要素を比較したらループを終了
        break
        
      } else if(vec[j] > val) {
        
        # j番目の値を1つ後に移動
        vec[j+1] <- vec[j]
        
      } else {
        
        # i番目の値以下ならループを終了
        break
      }
    }
    
    # i番目の値をj番目に挿入
    vec[j+1] <- val
  }
  
  # 数列を出力
  return(vec)
}

 数列を vec、数列の要素数を N とします。
 vec1 番目から N 番目まで(を i として)の要素を順番に入れ替えていきます。ただしアルゴリズム上、1番目の要素は入れ替えが行われないので、処理を省略します。
 veci 番目の要素 vec[i] の入れ替えでは、i 番目から 1 番目まで(を j として)の要素を順番に大小比較します。ただしアルゴリズム上、i番目の要素は入れ替えが行われないので、処理を省略します。

 i 番目の要素 vec[i]val としておき、valj 番目の要素 vec[j] を比較します。
 j 番目の値が大きければ、j 番目の値を j+1 番目の値に複製します。この操作が、j番目の要素を1つ後に移動するのに対応します。ただしこの時点では、j番目とj+1番目が同じ値です。
 j 番目の値が小さければ、break でループ処理を終了して、j+1 番目の要素を val で上書きします。この操作が、値の挿入に対応します。j+1番目の値は、1つ前の試行においてj+2番目に移動(複製)されています。
 全ての要素が大きかった場合は、j0 にしておき、j+1 (1)番目の要素を val で上書きします。

 アルゴリズムの詳しい解説は、本のcode 12.1などを参照してください。

 実装した関数を試してみます。

 要素数を指定して、数列を作成します。

# 要素数を指定
N <- 50

# 数列を生成
a <- sample(x = 1:(2*N), size = N, replace = TRUE)
a
##  [1] 98 69 41 52 52 68 71 66 85 16 54 24 87 89 98 74 68 17 17 56 67 86 36  7 67
## [26] 84 45 47 44 73 96 14 34 84 48 17  2 88 77 88 71 26  3 31 58 96 33 45 70  5

 乱数生成関数 sample() などで数値ベクトルを生成します。

 挿入ソートにより並べ替えます。

# ソート
insertion_sort(a)
##  [1]  2  3  5  7 14 16 17 17 17 24 26 31 33 34 36 41 44 45 45 47 48 52 52 54 56
## [26] 58 66 67 67 68 68 69 70 71 71 73 74 77 84 84 85 86 87 88 88 89 96 96 98 98

 目で確認するのはアレなので、組み込み関数のソート結果と比較して確認します。

# 確認
sum(!(insertion_sort(a) == sort(a)))
## [1] 0

 !TRUE, FALSE を反転させて sum() で合計することで、一致しない要素数を得られます。結果が 0 であれば全ての要素が一致したことが分かります。

ソートの可視化

 次は、挿入ソートのアルゴリズムをグラフで確認します。

 数列の初期値を作成して、グラフで確認します。

 要素数を指定して、数列を作成します。

# 要素数を指定
N <- 20

# 数列を生成
a <- rnorm(n = N, mean = 0, sd = 1) |> 
  round(digits = 1)
table(a)
## a
## -1.9 -1.5 -0.9 -0.6 -0.4 -0.2 -0.1  0.1  0.2  0.5  0.6  0.8  1.1  1.2 
##    1    1    1    1    2    1    2    2    3    2    1    1    1    1

 この例では、(負の値を含めるため)正規乱数生成関数 rnorm() で値を生成して、(重複を含めるため)丸め込み関数 round() で小数第一位にしておきます。

 数列をデータフレームに格納します。

# 要素数を取得
N <- length(a)

# 数列を格納
tmp_df <- tibble::tibble(
  iteration = 0,   # 試行回数
  id        = 1:N, # 元のインデックス
  index     = 1:N, # 各試行のインデックス
  value     = a    # 要素
)
tmp_df
## # A tibble: 20 × 4
##    iteration    id index value
##        <dbl> <int> <int> <dbl>
##  1         0     1     1   0.5
##  2         0     2     2   1.1
##  3         0     3     3  -0.6
##  4         0     4     4   1.2
##  5         0     5     5  -0.4
##  6         0     6     6   0.6
##  7         0     7     7  -1.5
##  8         0     8     8   0.2
##  9         0     9     9   0.5
## 10         0    10    10   0.1
## 11         0    11    11  -0.9
## 12         0    12    12   0.8
## 13         0    13    13  -1.9
## 14         0    14    14  -0.1
## 15         0    15    15  -0.2
## 16         0    16    16  -0.4
## 17         0    17    17  -0.1
## 18         0    18    18   0.2
## 19         0    19    19   0.2
## 20         0    20    20   0.1

 作図用に、インデックスなどの値とあわせて数列を格納します。数列のみが与えられている場合は、要素数を N とします。

 数列を棒グラフで確認します。

# 数列を作図
ggplot() + 
  geom_bar(data = tmp_df, 
           mapping = aes(x = index, y = value, fill = factor(value)), stat = "identity") + 
  theme(panel.grid.minor.x = element_blank()) + # 図の体裁
  labs(title = "numerical sequence", 
       subtitle = paste0("iteration: ", unique(tmp_df[["iteration"]])), 
       fill = "value", 
       x = "index", y = "value")

 この要素(バー)を昇順に並べ替えます。

 アルゴリズムに従いソートを行って、作図用のデータを作成します。

 1データずつ入れ替え処理を行い、数列を記録します。

# 作図用のオブジェクトを初期化
id_vec   <- 1:N
trace_df <- tmp_df

# 挿入ソート
for(i in 1:N) {
  if(i > 1) { # (初回は不要)
    
    # i番目の値を取得
    v <- a[i]
    
    # i番目から逆順に挿入位置jを探索:(j ≤ i)
    for(j in (i-1):0) { # (i番目は不要・0番目は挿入処理との兼ね合い用)
      
      # j番目とi番目の値を大小比較
      if(j <= 0) {
        
        # 全要素を比較したらループを終了
        break
        
      } else if(a[j] > v) {
        
        # j番目の値を1つ後に移動
        a[j+1] <- a[j]
        id_vec[1:N] <- replace(x = id_vec, list = j+1, values = id_vec[j])
        
      } else {
        
        # i番目の値以下ならループを終了
        break
      }
    }
    
    # i番目の値をj番目に挿入
    a[j+1] <- v
    id_vec[1:N] <- replace(x = id_vec, list = j+1, values = i)
  }
  
  # 数列を格納
  tmp_df <- tibble::tibble(
    iteration = i, 
    id        = id_vec, 
    index     = 1:N, 
    value     = a
  )
  
  # 数列を記録
  trace_df <- dplyr::bind_rows(trace_df, tmp_df)
  
  # 途中経過を表示
  #print(paste0("--- iteration: ", i, " ---"))
  #print(a)
}

 「ソートの実装」のときと同様にして、要素を入れ替えていきます。
 数列 a の初期値のインデックス(1 から N の整数)を id_vec としておき、a の要素の入れ替えに対応するように id_vec の要素も入れ替えます。
 試行ごとに、数列やインデックスをデータフレームに保存します。

 初期値のときのコードで、ソートした数列をグラフで確認します。

 昇順に並んでいるのが分かります。

 各試行の数列のグラフを並べて推移を確認します。

 挿入データの描画用のデータフレームを作成します。

# 挿入データを作成
target_df <- trace_df |> 
  dplyr::filter(iteration == id) # i番目の要素を抽出
target_df
## # A tibble: 20 × 4
##    iteration    id index value
##        <dbl> <int> <int> <dbl>
##  1         1     1     1   0.5
##  2         2     2     2   1.1
##  3         3     3     1  -0.6
##  4         4     4     4   1.2
##  5         5     5     2  -0.4
##  6         6     6     4   0.6
##  7         7     7     1  -1.5
##  8         8     8     4   0.2
##  9         9     9     6   0.5
## 10        10    10     4   0.1
## 11        11    11     2  -0.9
## 12        12    12    10   0.8
## 13        13    13     1  -1.9
## 14        14    14     6  -0.1
## 15        15    15     6  -0.2
## 16        16    16     6  -0.4
## 17        17    17     9  -0.1
## 18        18    18    12   0.2
## 19        19    19    13   0.2
## 20        20    20    11   0.1

 試行回数(iteration 列)と元のインデックス(id 列)が一致するデータ(行)を抽出します。

 試行ごとに数列のグラフを作成します。

# 全試行の数列を作図
ggplot() + 
  geom_bar(data = trace_df, 
           mapping = aes(x = index, y = value, fill = factor(value)), stat = "identity") + # 全ての要素
  geom_bar(data = target_df,
           mapping = aes(x = id, y = value), stat = "identity",
           color = "red", alpha = 0, linewidth = 0.5, linetype = "dashed") + # 挿入前のi番目の要素
  geom_bar(data = target_df,
           mapping = aes(x = index, y = value), stat = "identity",
           color = "red", alpha = 0, linewidth = 0.5, linetype = "dotted") + # 挿入後のi番目の要素
  geom_segment(data = target_df, 
               mapping = aes(x = id, y = value, xend = -Inf, yend = value, color = factor(value)), 
               linewidth = 0.5, linetype = "dashed", show.legend = FALSE) + # 挿入データまでの上限値
  facet_wrap(iteration ~ ., scales = "free_x", labeller = label_both) + # 試行ごとに分割
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "right") + 
  labs(title = "insertion sort", 
       fill = "value", 
       x = "index", y = "value")

 facet_wrap() で、試行(iteration 列)ごとに分割して棒グラフを描画します。
 挿入前の要素(バー)の位置は元のインデックス(id 列)、挿入後の位置は iteration 回目時点のインデックス(index 列)を使います。

 初期値を0回目として、N回目までの「入れ替え後の数列」をそれぞれ棒グラフ(バー)で表します。
 i回目の試行(i+1枚目のグラフ)において、挿入した要素を赤色の枠線のバー、またその値を水平線で示します。破線のバーが挿入前の位置、点線のバーが挿入後の位置です。
 挿入前の位置より左(小さいインデックス)の要素(バー)はソート済み、右(大きいインデックス)の要素(バー)は未処理のデータです。
 挿入後の位置より左の要素は挿入した要素よりも値(高さ)が小さく、右の要素は値が大きいのが分かります。

 各試行の数列のグラフを切り替えて推移を確認します。

 同様に、挿入データの描画用のデータフレームを作成します。

# 挿入データを作成
target_df <- trace_df |> 
    dplyr::filter((iteration + 1) == index) # i番目の要素を抽出
target_df
## # A tibble: 20 × 4
##    iteration    id index value
##        <dbl> <int> <int> <dbl>
##  1         0     1     1   0.5
##  2         1     2     2   1.1
##  3         2     3     3  -0.6
##  4         3     4     4   1.2
##  5         4     5     5  -0.4
##  6         5     6     6   0.6
##  7         6     7     7  -1.5
##  8         7     8     8   0.2
##  9         8     9     9   0.5
## 10         9    10    10   0.1
## 11        10    11    11  -0.9
## 12        11    12    12   0.8
## 13        12    13    13  -1.9
## 14        13    14    14  -0.1
## 15        14    15    15  -0.2
## 16        15    16    16  -0.4
## 17        16    17    17  -0.1
## 18        17    18    18   0.2
## 19        18    19    19   0.2
## 20        19    20    20   0.1

 こちらは、試行回数(iteration 列)と元のインデックス(id 列)が一致する要素の1試行前のデータを抽出します。

 重複ラベルの描画用のデータフレームを作成します。

# 重複ラベルを作成
dup_label_df <- trace_df |> 
  dplyr::arrange(iteration, id) |> # IDの割当用
  dplyr::group_by(iteration, value) |> # IDの割当用
  dplyr::mutate(
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, 
      true = paste0("(", dup_id, ")"), 
      false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup() |> 
  dplyr::arrange(iteration, index)
dup_label_df
## # A tibble: 420 × 7
##    iteration    id index value dup_id dup_num dup_label
##        <dbl> <int> <int> <dbl>  <int>   <int> <chr>    
##  1         0     1     1   0.5      1       2 "(1)"    
##  2         0     2     2   1.1      1       1 ""       
##  3         0     3     3  -0.6      1       1 ""       
##  4         0     4     4   1.2      1       1 ""       
##  5         0     5     5  -0.4      1       2 "(1)"    
##  6         0     6     6   0.6      1       1 ""       
##  7         0     7     7  -1.5      1       1 ""       
##  8         0     8     8   0.2      1       3 "(1)"    
##  9         0     9     9   0.5      2       2 "(2)"    
## 10         0    10    10   0.1      1       2 "(1)"    
## # ℹ 410 more rows

 重複している値の要素にのみ、それぞれ通し番号のラベルを作成します。

 推移のアニメーションを作成します。

# ソートのアニメーションを作図
graph <- ggplot() + 
  geom_bar(data = trace_df, 
           mapping = aes(x = index, y = value, fill = factor(value), group = factor(id)), stat = "identity") + # 全要素
  geom_bar(data = target_df,
           mapping = aes(x = index, y = value, group = factor(id)), stat = "identity",
           color = "red", alpha = 0, linewidth = 1, linetype = "dashed") + # 挿入データ
  geom_segment(data = target_df,
               mapping = aes(x = index, y = value, xend = -2, yend = value,
                             color = factor(value), group = factor(id)),
               linewidth = 1, linetype = "dashed", show.legend = FALSE) + # 挿入データまでの上限値
  geom_text(data = trace_df, 
            mapping = aes(x = index, y = 0, label = as.character(value), group = factor(id)), 
            vjust = -0.5, size = 5) + # 要素ラベル
  geom_text(data = dup_label_df, 
            mapping = aes(x = index, y = 0, label = dup_label, group = factor(id)), 
            vjust = 1, size = 4) + # 重複ラベル
  gganimate::transition_states(states = iteration, transition_length = 9, state_length = 1, wrap = FALSE) + # フレーム遷移
  gganimate::ease_aes("cubic-in-out") + # 遷移の緩急
  coord_cartesian(xlim = c(0, N+1)) + # (上限値の線用)
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "insertion sort", 
       subtitle = "iteration: {next_state}", 
       fill = "value", 
       x = "index", y = "value")

# 1試行当たりのフレーム数を指定
s <- 20

# gif画像を作成
gganimate::animate(
  plot = graph, 
  nframes = (N+1 + 2)*s, start_pause = s, end_pause = s, fps = 20, 
  width = 800, height = 600, 
  renderer = gganimate::gifski_renderer()
)

 transition_states() のフレーム制御の引数 states に試行回数列 iteration を指定して、これまでと同様に作図します。
 animate()plot 引数にグラフオブジェクト、nframes 引数にフレーム数を指定して、gif画像を作成します。
 初期値を含むため試行回数は N+1 です。transition_states() でアニメーションを作成すると、フレーム間の値を補完して、状態の遷移を描画します。遷移のフレームを s として、試行回数の s 倍の値を nframes 引数に指定します。

 左(1番目)の要素から順番にソートされるのを確認できます。また、重複要素の順序が入れ替わらない(安定性・stable性)も確認できます。

 ここでの試行回数(iterationの値)は、この可視化方法における試行番号(作図処理上の値)であり他の手法と対応していません。

 この記事では、挿入ソートを実装しました。次の記事では、マージソートを実装します。

12.4 マージソートの実装と可視化

 マージソート(merge sort)を実装します。また、ソートの推移のアニメーションを作成して、アルゴリズムを確認します。

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(gganimate)

 この記事では、基本的に パッケージ名::関数名() の記法を使っているので、パッケージを読み込む必要はありません。ただし、作図コードについてはパッケージ名を省略しているので、ggplot2 を読み込む必要があります。
 また、ネイティブパイプ演算子 |> を使っています。magrittr パッケージのパイプ演算子 %>% に置き換えても処理できますが、その場合は magrittr を読み込む必要があります。

 「ソートの可視化」において、gganimate パッケージを利用してアニメーション(gif画像)を作成します。ソートの実装自体には使いません。

ソートの実装

 まずは、マージソートのアルゴリズムを関数として実装します。

 マージソートを実装します。

# マージソートの実装
merge_sort <- function(vec) {
  
  # 入替範囲の要素数を取得
  l <- 1
  r <- length(vec)
  
  # 要素数が1なら再帰処理を終了
  if(r == 1) return(vec)
  
  # 分割位置mを計算:(l ≤ m ≤ r)
  m <- l + (r - l) %/% 2
  
  # 分割してそれぞれソート
  buf <- c(
    merge_sort(vec[l:m]),
    merge_sort(vec[(m+1):r]) |>
      rev()
  )
  
  # 入替範囲内のインデックスを初期化
  idx_l <- 1
  idx_r <- 0
  
  # 要素ごとに処理
  for(i in 1:r) {
    
    # 小さい方の要素を取り出して併合
    if(buf[idx_l] <= buf[r-idx_r]) { # 左側が小さい場合
      
      # 左側の要素をj番目に格納
      vec[i] <- buf[idx_l]
      
      # インデックスを更新
      idx_l  <- idx_l + 1
      
    } else { # 右側が小さい場合
      
      # 右側の要素をj番目に格納
      vec[i] <- buf[r-idx_r]
      
      # インデックスを更新
      idx_r  <- idx_r + 1
    }
  }
  
  # 数列を出力
  return(vec)
}

 ソート対象の数列のうち、分割された一部の要素を vec として受け取ります。
 vec の要素数を r とします。また、l1 としておきます。(数列全体を見た際のlからrの範囲に対応した変数名です。ここでは一部だけを扱っているので1とnの方が分かりやすいかもしれません。)
 vec (l から r の範囲)をさらに分割する位置(インデックス)を m とします。

 「l 番目から m 番目までの要素」と「m+1 番目から r 番目までの要素」に分割してそれぞれ merge_sort() でソートします。
 2つに分割したうち、「前(1から)の方はそのまま」で「後(m+1から)の方は逆順」にして結合したベクトルを一時変数 buf とします。1 から m 番目までは昇順、m+1 から r 番目までは降順に並んだ状態です。
 頭と尻の要素の大小を(勝ち抜け戦の要領で)比較して、小さい順に取り出して vec に格納していきます。

 merge_sort() を実行すると内部でさらに(全ての分割範囲が1になるまで再帰的に) merge_sort() が実行され、分割と併合(ソート)が繰り返されて、全体がソートされます。

 アルゴリズムの詳しい解説は、本のcode 12.2などを参照してください。

 実装した関数を試してみます。

 要素数を指定して、数列を作成します。

# 要素数を指定
N <- 50

# 数列を生成
a <- sample(x = 1:(2*N), size = N, replace = TRUE)
a
##  [1] 55 34 85 37 96 34 77 36 96  5  5 29  2 73 40 14 32 20 93 58 87 54 14  5 38
## [26] 87 35 43 31 75 81 78 95 51 63  7 75 78 49 15 83 53 10 80 13 62 17 70 37 80

 乱数生成関数 sample() などで数値ベクトルを生成します。

 マージソートにより並べ替えます。

# ソート
merge_sort(a)
##  [1]  2  5  5  5  7 10 13 14 14 15 17 20 29 31 32 34 34 35 36 37 37 38 40 43 49
## [26] 51 53 54 55 58 62 63 70 73 75 75 77 78 78 80 80 81 83 85 87 87 93 95 96 96

 目で確認するのはアレなので、組み込み関数のソート結果と比較して確認します。

# 確認
sum(!(merge_sort(a) == sort(a)))
## [1] 0

 !TRUE, FALSE を反転させて sum() で合計することで、一致しない要素数を得られます。結果が 0 であれば全ての要素が一致したことが分かります。

ソートの可視化

 次は、マージソートのアルゴリズムをグラフで確認します。

 数列の初期値を作成して、グラフで確認します。

 要素数を指定して、数列を作成します。

# 要素数を指定
N <- 35

# 数列を生成
a <- rnorm(n = N, mean = 0, sd = 1) |> 
  round(digits = 1)
table(a)
## a
## -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.1  0.1  0.3  0.4 
##    1    2    1    1    1    1    1    2    3    1    4    3    1    1    2    1 
##  0.5  0.8    1  1.2  1.4  1.7  1.8 
##    2    1    1    2    1    1    1

 この例では、(負の値を含めるため)正規乱数生成関数 rnorm() で値を生成して、(重複を含めるため)丸め込み関数 round() で小数第一位にしておきます。

 数列をデータフレームに格納します。

# 要素数を取得
N <- length(a)

# 数列を格納
tmp_df <- tibble::tibble(
  iteration = 0,   # 試行回数
  id        = 1:N, # 元のインデックス
  index     = 1:N, # 各試行のインデックス
  value     = a    # 要素
)
tmp_df
## # A tibble: 35 × 4
##    iteration    id index value
##        <dbl> <int> <int> <dbl>
##  1         0     1     1  -0.7
##  2         0     2     2  -0.8
##  3         0     3     3  -0.1
##  4         0     4     4   0.8
##  5         0     5     5   0.1
##  6         0     6     6  -0.6
##  7         0     7     7   1.2
##  8         0     8     8   0.3
##  9         0     9     9  -0.7
## 10         0    10    10  -1.1
## # ℹ 25 more rows

 作図用に、インデックスなどの値とあわせて数列を格納します。数列のみが与えられている場合は、要素数を N とします。

 数列を棒グラフで確認します。

# 数列を作図
ggplot() + 
  geom_bar(data = tmp_df, 
           mapping = aes(x = index, y = value, fill = factor(value)), stat = "identity") + 
  theme(panel.grid.minor.x = element_blank()) + 
  labs(title = "numerical sequence", 
       subtitle = paste0("iteration: ", unique(tmp_df[["iteration"]])), 
       fill = "value", 
       x = "index", y = "value")

 この要素(バー)を昇順に並べ替えます。

 アルゴリズムに従いソートを行って、作図用のデータを作成します。
 この例では、数列の全ての範囲に対して分割または併合を行う処理を1試行とします。再帰的な処理を行うと次の試行になります。ただし実際(関数内部)の処理では、1試行の処理が全て終わる前に再帰処理が行われることもあります。

 マージソートでは、1要素ずつになるまで数列を二分割していき、分割した数列の一部を(分割時とは逆方向に)併合しながら昇順に並べ替えていくのでした。つまり、要素の入れ替えを行う前に数列を分割しておく必要があります。
 そこで、数列の分割に対応するインデックスを作成する関数を定義しておきます。

# インデックスの分割の実装
trace_idx <- function(mat = NA, iter = 0, left = 1, right = 1) {
  
  # 分割範囲を取得
  i <- iter
  l <- left
  r <- right
  
  # 試行回数をカウント
  if(i == 0) {
    
    # 初回ならマトリクスを初期化
    mat <- matrix(data = NA, nrow = 1, ncol = r)
    i <- 1
    
  } else if(l == 1) {
    
    # 左端の範囲なら行を追加
    mat <- rbind(
      mat, 
      matrix(data = NA, nrow = 1, ncol = ncol(mat))
    )
    i <- i + 1
    
  } else {
    
    i <- i + 1
  }
  
  # 分割位置mを計算:(l ≤ m ≤ r)
  m <- l + (r - l) %/% 2
  
  # 分割範囲のインデックスを格納
  mat[i, l:r] <- 1:(r-l+1)
  
  # 範囲が1なら再帰処理を終了
  if(l == r) return(mat)
  
  # インデックスを取得
  mat <- trace_idx(mat, iter = i, left = l, right = m)
  mat <- trace_idx(mat, iter = i, left = (m+1), right = r)
  
  # インデックスを出力
  return(mat)
}

 (分割の仕方自体は重要ではないと思いますが、)merge_sort() の内部と同じ挙動になるように(頑張って)処理します。
 試行ごとにマトリクスに行を追加します。l1 のとき、つまりソート対象の数列の左端の要素を含む範囲のとき、新たな試行とみなします。
 分割範囲ごとに通し番号のインデックスを割り当てますが、使うのは 1 の位置(列番号)のみです。

 実装した関数を使って、分割した際のインデックスを作成します。

# 分割インデックスを作成
trace_idx_mat <- trace_idx(mat = NA, iter = 0, left = 1, right = N)
trace_idx_mat[, 1:5]
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    2    3    4    5
## [2,]    1    2    3    4    5
## [3,]    1    2    3    4    5
## [4,]    1    2    3    4    5
## [5,]    1    2    3    1    2
## [6,]    1    2    1    1    1
## [7,]    1    1   NA   NA   NA

 分割された範囲の要素数が1になると、その範囲では再帰処理が行われません。そのため、他の範囲で再帰処理が行われると、既に要素数が1の範囲の次の行は欠損値 NA になります。
 trace_idx() が出力したマトリクスに対して、is.na()TRUE になる要素(欠損値)に 1 を代入します。

# 欠損値を補完
trace_idx_mat[is.na(trace_idx_mat)] <- 1
trace_idx_mat[, 1:5]
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    2    3    4    5
## [2,]    1    2    3    4    5
## [3,]    1    2    3    4    5
## [4,]    1    2    3    4    5
## [5,]    1    2    3    1    2
## [6,]    1    2    1    1    1
## [7,]    1    1    1    1    1

 数列の要素数に応じて、分割(入替)の範囲が決まります。
 最後の行は、全ての要素が 1 になり、1要素ずつに分割された状態を表し、また併合されていない状態(数列の初期値)を表します。

 1試行ずつ入れ替え処理を行い、数列を記録します。

# 試行回数を取得:(初期値は除く)
max_iter <- nrow(trace_idx_mat) - 1

# 作図用のオブジェクトを初期化
id_vec   <- 1:N
trace_df <- tmp_df

# マージソート
for(i in 1:max_iter) {
  
  # 分割範囲のインデックスを作成
  l_idx_vec <- which(trace_idx_mat[max_iter-i+1, ] == 1)
  r_idx_vec <- c(l_idx_vec[-1]-1, N)
  m_idx_vec <- l_idx_vec + (r_idx_vec - l_idx_vec) %/% 2
  
  # 分割数を取得
  max_k <- length(l_idx_vec)
  
  # 分割範囲ごとに処理
  for(k in 1:max_k) {
    
    # 分割範囲を取得
    l <- l_idx_vec[k]
    r <- r_idx_vec[k]
    m <- m_idx_vec[k]
    
    # 前回の要素IDを保存
    old_id_vec <- id_vec
    
    # 数列を分割
    if(l == r) { # 要素数が1の場合
      buf <- a[l]
    } else {
      buf <- c(a[l:m], rev(a[(m+1):r]))
    }
    
    # 分割範囲内のインデックスを初期化
    idx_l <- 1
    idx_r <- 0
    max_r <- length(buf)
    
    # 要素ごとに処理
    for(j in l:r) {
      
      # 小さい方の要素を取り出して併合
      if(buf[idx_l] <= buf[max_r-idx_r]) {
        
        # 左側の要素をj番目に格納
        a[j] <- buf[idx_l]
        id_vec[1:N] <- replace(x = id_vec, list = j, values = old_id_vec[l+idx_l-1])
        
        # インデックスを更新
        idx_l <- idx_l + 1
        
      } else {
        
        # 右側の要素をj番目に格納
        a[j] <- buf[max_r-idx_r]
        id_vec[1:N] <- replace(x = id_vec, list = j, values = old_id_vec[m+idx_r+1])
        
        # インデックスを更新
        idx_r <- idx_r + 1
      }
    }
  }
  
  # 数列を格納
  tmp_df <- tibble::tibble(
    iteration = i, 
    id        = id_vec, 
    index     = 1:N, 
    value     = a
  )
  
  # 数列を記録
  trace_df <- dplyr::bind_rows(trace_df, tmp_df)

  # 途中経過を表示
  #print(paste0("--- iteration: ", i, " ---"))
  #print(a)
}

 分割した範囲内のインデックス trace_idx_mat の最後の行(初期値に対応するインデックス)を除き、後の行から順番に取り出して処理します。そのため、trace_idx_mat の行数-1を試行回数の最大値 max_iter とします。
 trace_idx_mat の各行について、値が 1 の要素の列インデックスが、分割範囲の左端のインデックス l_idx_vec に対応します。
 l_idx_vec の最初の要素(1)を除き、最後に要素数(N)を追加したベクトルが、分割範囲の右端のインデックス r_idx_vec に対応します。
 l_idx_vec, r_idx_vec から分割点のインデックス m_idx_vec を作成します。m_idx_vec + 1l_idx_vec の値が併合における1試行前(分割だと1試行後)の l_idx_vec の値です。
 *_idx_vec の要素数が各試行の分割数に対応します。

 分割範囲ごとに分割範囲と分割位置のインデックスと取り出して、「ソートの実装」のときと同様にして、要素を入れ替えていきます。
 数列 a の初期値のインデックス(1 から N の整数)を id_vec としておき、a の要素の入れ替えに対応するように id_vec の要素も入れ替えます。

 試行ごとに、数列やインデックスをデータフレームに保存します。

 初期値のときのコードで、ソートした数列をグラフで確認します。

 昇順に並んでいるのが分かります。

 各試行の数列のグラフを並べて推移を確認します。

 分割範囲(入替範囲)の描画用のデータフレームを作成します。

# 入替範囲を作成
range_df <- trace_idx_mat[(max_iter+1):1, ] |> # (逆順に並べ替え)
  tibble::as_tibble(.name_repair = NULL) |> 
  tibble::add_column(iteration = 0:max_iter) |> # 試行回数列を追加
  tidyr::pivot_longer(
    cols = !iteration, 
    names_to = "index", 
    names_prefix = "V", 
    names_transform = list(index = as.numeric), 
    values_to = "split_flag"
  ) |> # 分割インデックス列をまとめる
  dplyr::filter(split_flag == 1) |> # 分割位置を抽出
  dplyr::arrange(iteration, index) |> # 分割範囲の作成用
  dplyr::group_by(iteration) |> # 分割範囲の作成用
  dplyr::mutate(
    # 分割範囲を作成
    left  = index, 
    right = dplyr::lead(index, n = 1, default = N+1), 
    # タイル用
    x      = 0.5 * (right + left - 1), 
    y      = 0.5 * (max(a) + min(c(0, a))), 
    width  = right - left, 
    height = max(a) - min(c(0, a)) + 1
  ) |> 
  dplyr::ungroup()
range_df
## # A tibble: 98 × 9
##    iteration index split_flag  left right     x     y width height
##        <int> <dbl>      <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
##  1         0     1          1     1     2     1   0.1     1    4.4
##  2         0     2          1     2     3     2   0.1     1    4.4
##  3         0     3          1     3     4     3   0.1     1    4.4
##  4         0     4          1     4     5     4   0.1     1    4.4
##  5         0     5          1     5     6     5   0.1     1    4.4
##  6         0     6          1     6     7     6   0.1     1    4.4
##  7         0     7          1     7     8     7   0.1     1    4.4
##  8         0     8          1     8     9     8   0.1     1    4.4
##  9         0     9          1     9    10     9   0.1     1    4.4
## 10         0    10          1    10    11    10   0.1     1    4.4
## # ℹ 88 more rows

 分割した範囲内のインデックス trace_idx_mat をデータフレームに変換して、pivot_longer() で全ての値を1つの列にまとめます。列名(を加工した値)は数列全体のインデックスに対応します。
 trace_idx_mat の行番号の逆順に 0 から max_iter までの整数を割り当てて、試行回数列とします。

 分割範囲内インデックスが 1 の要素(行)を取り出します。抽出された要素の全体のインデックスが、分割範囲の左端の位置lに対応します。また、1 を除き N+1 を追加したインデックスが右端の位置rに対応します。
 そこで、lead() で1行前にズラして、最後に(default 引数で) N+1 を追加した列を作成します。

 試行ごとに数列のグラフを作成します。

# 全試行の数列を作図
ggplot() + 
  geom_bar(data = trace_df, 
           mapping = aes(x = index, y = value, fill = factor(value)), stat = "identity") + # 全ての要素
  geom_tile(data = range_df,
            mapping = aes(x = x, y = y, width = width, height = height),
            color = "red", alpha = 0, linewidth = 0.6, linetype = "dashed") + # 入替範囲
  facet_wrap(iteration ~ ., scales = "free_x", labeller = label_both) + # 試行ごとに分割
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "merge sort", 
       fill = "value", 
       x = "index", y = "value")

 facet_wrap() で、試行(iteration 列)ごとに分割して棒グラフを描画します。
 要素の入れ替えが行われた(分割された)範囲を geom_tile() で描画します。(geom_bar() を使って範囲を描画しようとすると、描画自体はできるのですが、width 引数について警告文が出たので避けました。原因は不明です。)

 初期値を0回目として、各試行における「入れ替え後の数列」をそれぞれ棒グラフ(バー)で表します。
 i回目の試行(i+1枚目のグラフ)において、入れ替えの行われた(分割された)範囲を赤色の破線の枠で示します。

 各グラフを見ると、枠内の要素(バー)が昇順に並んでいるのが分かります。また次のグラフを見ると、隣り合う枠(範囲)を併合しながら全体がソートされるのが分かります。

 各試行の数列のグラフを切り替えて推移を確認します。

 同様に、入替範囲の描画用のデータフレームを作成します。

# 入替範囲を作成
range_df <- trace_idx_mat[max_iter:1, ] |> # (最後を除き逆順に並べ替え)
  tibble::as_tibble(.name_repair = NULL) |> 
  tibble::add_column(iteration = 0:(max_iter-1)) |> # 試行回数列を追加
  tidyr::pivot_longer(
    cols = !iteration, 
    names_to = "index", 
    names_prefix = "V", 
    names_transform = list(index = as.numeric), 
    values_to = "split_flag"
  ) |> # 分割インデックス列をまとめる
  dplyr::filter(split_flag == 1) |> # 分割位置を抽出
  dplyr::arrange(iteration, index) |> # 分割範囲の作成用
  dplyr::group_by(iteration) |> # 分割範囲の作成用
  dplyr::mutate(
    # 分割範囲を作成
    left  = index, 
    right = dplyr::lead(index, n = 1, default = N+1), 
    # タイル用
    x      = 0.5 * (right + left - 1), 
    y      = 0.5 * (max(a) + min(c(0, a))), 
    width  = right - left, 
    height = max(a) - min(c(0, a)) + 1
  ) |> 
  dplyr::ungroup() |> 
  dplyr::mutate(
    id = dplyr::row_number() # 通し番号
  )
range_df
## # A tibble: 63 × 10
##    iteration index split_flag  left right     x     y width height    id
##        <int> <dbl>      <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <int>
##  1         0     1          1     1     3   1.5   0.1     2    4.4     1
##  2         0     3          1     3     4   3     0.1     1    4.4     2
##  3         0     4          1     4     5   4     0.1     1    4.4     3
##  4         0     5          1     5     6   5     0.1     1    4.4     4
##  5         0     6          1     6     7   6     0.1     1    4.4     5
##  6         0     7          1     7     8   7     0.1     1    4.4     6
##  7         0     8          1     8     9   8     0.1     1    4.4     7
##  8         0     9          1     9    10   9     0.1     1    4.4     8
##  9         0    10          1    10    12  10.5   0.1     2    4.4     9
## 10         0    12          1    12    13  12     0.1     1    4.4    10
## # ℹ 53 more rows

 こちらは、1試行分タイミングを早めてデータを作成します。

 重複ラベルの描画用のデータフレームを作成します。

# 重複ラベルを作成
dup_label_df <- trace_df |> 
  dplyr::arrange(iteration, id) |> # IDの割当用
  dplyr::group_by(iteration, value) |> # IDの割当用
  dplyr::mutate(
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, 
      true = paste0("(", dup_id, ")"), 
      false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup() |> 
  dplyr::arrange(iteration, index)
dup_label_df
## # A tibble: 245 × 7
##    iteration    id index value dup_id dup_num dup_label
##        <dbl> <int> <int> <dbl>  <int>   <int> <chr>    
##  1         0     1     1  -0.7      1       2 "(1)"    
##  2         0     2     2  -0.8      1       1 ""       
##  3         0     3     3  -0.1      1       1 ""       
##  4         0     4     4   0.8      1       1 ""       
##  5         0     5     5   0.1      1       1 ""       
##  6         0     6     6  -0.6      1       3 "(1)"    
##  7         0     7     7   1.2      1       2 "(1)"    
##  8         0     8     8   0.3      1       2 "(1)"    
##  9         0     9     9  -0.7      2       2 "(2)"    
## 10         0    10    10  -1.1      1       1 ""       
## # ℹ 235 more rows

 重複している値の要素にのみ、それぞれ通し番号のラベルを作成します。

 推移のアニメーションを作成します。

# ソートのアニメーションを作図
graph <- ggplot() + 
  geom_bar(data = trace_df, 
           mapping = aes(x = index, y = value, fill = factor(value)), stat = "identity") + # 全ての要素
  geom_tile(data = range_df,
            mapping = aes(x = x, y = y, width = width, height = height, group = factor(id)),
            color = "red", alpha = 0, linewidth = 1, linetype = "dashed") + # 入替範囲
  geom_text(data = trace_df, 
            mapping = aes(x = index, y = 0, label = as.character(value), group = factor(id)), 
            vjust = -0.5, size = 5) + # 要素ラベル
  geom_text(data = dup_label_df, 
            mapping = aes(x = index, y = 0, label = dup_label, group = factor(id)), 
            vjust = 1, size = 4) + # 重複ラベル
  gganimate::transition_states(states = iteration, transition_length = 9, state_length = 1, wrap = FALSE) + # フレーム遷移
  gganimate::ease_aes("cubic-in-out") + # 遷移の緩急
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "merge sort", 
       subtitle = "iteration: {next_state}", 
       fill = "value", 
       x = "index", y = "value")

# 遷移フレーム数を指定
s <- 20

# gif画像を作成
gganimate::animate(
  plot = graph, 
  nframes = (max_iter+1 + 2)*s, start_pause = s, end_pause = s, fps = 20, 
  width = 1200, height = 900, 
  renderer = gganimate::gifski_renderer()
)

 transition_states() のフレーム制御の引数 states に試行回数列 iteration を指定して、これまでと同様に作図します。
 animate()plot 引数にグラフオブジェクト、nframes 引数にフレーム数を指定して、gif画像を作成します。
 初期値を含むため試行回数は max_iter+1 (trace_idx_mat の行数)です。transition_states() でアニメーションを作成すると、フレーム間の値を補完して、状態の遷移を描画します。遷移のフレームを s として、試行回数の s 倍の値を nframes 引数に指定します。

 分割された範囲が併合する度にソートされるのを確認できます。
 また、重複要素の順序が入れ替わらない(安定性・stable性)のも分かります。

 ここでの試行回数(iterationの値)は、この可視化方法における試行番号(作図処理上の値)であり他の手法と対応していません。

 この記事では、マージソートを実装しました。次の記事では、クイックソートを実装します。

12.5 クイックソートの実装と可視化

 クイックソート(quick sort)を実装します。また、ソートの推移のアニメーションを作成して、アルゴリズムを確認します。

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(patchwork)
library(gganimate)

 この記事では、基本的に パッケージ名::関数名() の記法を使っているので、パッケージを読み込む必要はありません。ただし、作図コードについてはパッケージ名を省略しているので、ggplot2 を読み込む必要があります。
 また、ネイティブパイプ演算子 |> を使っています。magrittr パッケージのパイプ演算子 %>% に置き換えても処理できますが、その場合は magrittr を読み込む必要があります。

 「ソートの可視化」において、patchwork パッケージを利用して複数のグラフを並べて描画、gganimate パッケージを利用してアニメーション(gif画像)を作成します。ソートの実装自体には使いません。

ソートの実装

 まずは、クイックソートのアルゴリズムを関数として実装します。

 クイックソートを実装します。

# クイックソートの実装
quick_sort <- function(vec) {
  
  # 入替範囲の要素数を取得
  l <- 1
  r <- length(vec)
  
  # 要素数が1なら再帰処理を終了
  if(r == 1) return(vec)
  
  # ピボットを設定
  p_idx <- (l + r) %/% 2 # (中点)
  p_val <- vec[p_idx]
  
  # ピボットを最後尾と入替
  vec[p_idx] <- vec[r]
  vec[r]     <- p_val
  
  # ピボット未満の値を探索:(j < r)
  i <- l
  for(j in l:(r-1)) { # (ピボットを除く)
    if(vec[j] < p_val) {
      
      # i番目とj番目の値を入替
      vec[l:r] <- replace(x = vec, list = c(i, j), values = vec[c(j, i)])
      
      # 境界位置を更新
      i <- i + 1
    }
  }
  
  # ピボットを境界と入替
  vec[r] <- vec[i]
  vec[i] <- p_val
  
  # ピボット前後で分割してソート
  vec[l:(i-1)] <- quick_sort(vec[l:(i-1)])
  if(i < r) { # (境界が最後尾なら不要)
    vec[(i+1):r] <- quick_sort(vec[(i+1):r])
  }
  
  # 数列を出力
  return(vec)
}

 ソート対象の数列のうち、分割された一部(入替範囲内)の要素を vec として受け取ります。
 vec の要素数を r とします。また、l1 としておきます。(数列全体を見た際のlからrの範囲に対応した変数名です。ここでは一部だけを扱っているので1とnの方が分かりやすいかもしれません。)
 vec (l から r の範囲)の要素からピボットを1つ決めます。この例では中点とします。ピボットのインデックスを p_idx、値を p_val とします。

 範囲内の最後(r番目)の要素とピボットの要素を入れ替えて、範囲内の一番後ろにピボットの要素を配置します。
 範囲内の最初(l番目)の要素から順番にピボットの値と比較して、小さければ i 番目の要素と入れ換えて、範囲内の前方にピボット未満の要素を配置していきます。要素を入れ替える(ピボット未満の要素が見付かる)と、i1 を加えます。i の初期値を l にしておくことで、i は「前方に集められたピボット未満の要素」の次の要素のインデックスになります。
 全ての要素を比較したら、i 番目の要素とピボットの要素を入れ替えます。入れ替え後のピボット位置は、ピボット未満の要素とピボット以上の要素の境界になります。

 「l 番目から i-1 番目まで(境界より前)の要素」と「i+1 番目から r 番目まで(境界より後)の要素」に分割してそれぞれ quick_sort() でソートします。
 quick_sort() を実行すると内部でさらに(全ての分割範囲が1になるまで再帰的に) quick_sort() が実行され、ピボットの設定と分割(ソート)が繰り返されて、全体がソートされます。

 アルゴリズムの詳しい解説は、本のcode 12.3などを参照してください。

 実装した関数を試してみます。

 要素数を指定して、数列を作成します。

# 要素数を指定
N <- 50

# 数列を生成
a <- sample(x = 1:(2*N), size = N, replace = TRUE)
a
##  [1]  5 65 31 59 73 85 21 66 26 33 61 81 51 18 73 64 62 59 85 17  9 40 33  9 68
## [26] 69 21 71 31 78 56 96 30 81 82 61 23 52 93 18 73 20 12 21 90 35  5 78 12 55

 乱数生成関数 sample() などで数値ベクトルを生成します。

 クイックソートにより並べ替えます。

# ソート
quick_sort(a)
##  [1]  5  5  9  9 12 12 17 18 18 20 21 21 21 23 26 30 31 31 33 33 35 40 51 52 55
## [26] 56 59 59 61 61 62 64 65 66 68 69 71 73 73 73 78 78 81 81 82 85 85 90 93 96

 目で確認するのはアレなので、組み込み関数のソート結果と比較して確認します。

# 確認
sum(!(quick_sort(a) == sort(a)))
## [1] 0

 !TRUE, FALSE を反転させて sum() で合計することで、一致しない要素数を得られます。結果が 0 であれば全ての要素が一致したことが分かります。

ソートの可視化

 次は、クイックソートのアルゴリズムをグラフで確認します。

 数列の初期値を作成して、グラフで確認します。

 要素数を指定して、数列を作成します。

# 要素数を指定
N <- 50

# 数列を生成
a <- rnorm(n = N, mean = 0, sd = 2) |> 
  round(digits = 1)
table(a)
## a
## -3.8 -3.6 -3.3   -3 -2.9 -2.7 -2.3 -2.2 -1.6 -1.5 -0.9 -0.7 -0.5 -0.4 -0.3 -0.2 
##    1    1    1    1    1    1    1    3    1    1    2    1    1    1    1    3 
## -0.1    0  0.1  0.2  0.3  0.5  0.7  0.8    1  1.1  1.2  1.3  2.2  2.3  2.9    3 
##    3    1    1    3    1    2    2    1    1    3    1    1    1    1    1    1 
##  3.2  3.3  4.6  5.6 
##    1    2    1    1

 この例では、(負の値を含めるため)正規乱数生成関数 rnorm() で値を生成して、(重複を含めるため)丸め込み関数 round() で小数第一位にしておきます。

 数列をデータフレームに格納します。

# 要素数を取得
N <- length(a)

# 数列の初期値を格納
tmp_df <- tibble::tibble(
  iteration = 0,   # 試行回数
  id        = 1:N, # 元のインデックス
  index     = 1:N, # 各試行のインデックス
  value     = a,   # 数列
  pivot_flag       = FALSE, # 各試行のピボット
  trace_pivot_flag = FALSE  # 各試行までのピボット
)
tmp_df
## # A tibble: 50 × 6
##    iteration    id index value pivot_flag trace_pivot_flag
##        <dbl> <int> <int> <dbl> <lgl>      <lgl>           
##  1         0     1     1   0   FALSE      FALSE           
##  2         0     2     2   4.6 FALSE      FALSE           
##  3         0     3     3  -3   FALSE      FALSE           
##  4         0     4     4   0.5 FALSE      FALSE           
##  5         0     5     5  -2.2 FALSE      FALSE           
##  6         0     6     6   0.2 FALSE      FALSE           
##  7         0     7     7  -1.5 FALSE      FALSE           
##  8         0     8     8   1.2 FALSE      FALSE           
##  9         0     9     9   2.9 FALSE      FALSE           
## 10         0    10    10  -2.2 FALSE      FALSE           
## # ℹ 40 more rows

 作図用に、インデックスなどの値とあわせて数列を格納します。数列のみが与えられている場合は、要素数を N とします。

 数列を棒グラフで確認します。

# 数列を作図
ggplot() + 
  geom_bar(data = tmp_df, 
           mapping = aes(x = index, y = value, fill = factor(value)), stat = "identity") + 
  theme(panel.grid.minor.x = element_blank()) + 
  labs(title = "numerical sequence", 
       subtitle = paste0("iteration: ", unique(tmp_df[["iteration"]])), 
       fill = "value", 
       x = "index", y = "value")

 この要素(バー)を昇順に並べ替えます。

 アルゴリズムに従いソートを行って、作図用のデータを作成します。
 この例では、数列の全ての範囲に対してピボットの設定と分割を行う処理を1試行とします。再帰的な処理を行うと次の試行になります。ただし実際(関数内部)の処理では、1試行の処理が全て終わる前に再帰処理が行われることもあります。

 数列の入れ替え時に、値と共にインデックスも入れ替えて出力する関数を定義しておきます。

# ピボットの設定の実装
pivot_split <- function(value, index, left, right) {
  
  # 全体の要素数を取得
  val_vec <- value
  idx_vec <- index
  N       <- length(val_vec)
  
  # 分割範囲を取得
  l <- left
  r <- right
  
  # 分割範囲の要素数が1なら終了
  if(r - l < 1) {
    
    # リストに格納
    lt <- list(value = val_vec, index = idx_vec, new_index = l, pivot_index = l)
    
    # 結果を出力
    return(lt)
  }
  
  # ピボットを設定
  p_idx <- (l + r) %/% 2 # (中点で固定)
  p_val <- val_vec[p_idx]
  
  # ピボットを最後尾と入替
  val_vec[1:N] <- replace(x = val_vec, list = c(p_idx, r), values = c(val_vec[r], p_val))
  idx_vec[1:N] <- replace(x = idx_vec, list = c(p_idx, r), values = idx_vec[c(r, p_idx)])
  
  # ピボット未満の値を探索:(j < r)
  i <- l
  for(j in l:(r-1)) { # (ピボットを除く)
    if(val_vec[j] < p_val) {
      
      # j番目を境界と入替
      val_vec[1:N] <- replace(x = val_vec, list = c(i, j), values = val_vec[c(j, i)])
      idx_vec[1:N] <- replace(x = idx_vec, list = c(i, j), values = idx_vec[c(j, i)])
      
      # 境界位置をカウント
      i <- i + 1
    }
  }
  
  # ピボットを境界と入替
  val_vec[1:N] <- replace(x = val_vec, list = c(r, i), values = c(val_vec[i], p_val))
  idx_vec[1:N] <- replace(x = idx_vec, list = c(r, i), values = idx_vec[c(i, r)])
  
  # リストに格納
  lt <- list(value = val_vec, index = idx_vec, new_index = i, old_index = p_idx)
  
  # 結果を出力
  return(lt)
}

 数列全体と入替範囲のインデックスを受け取り、quick_sort() と同様にピボットの設定と入れ替えを行い、作図用の値を出力します。再帰処理は行いません。
 作図用の値として、「数列の値・インデックス」のベクトルと「ピボットの入れ替え前後のインデックス」の値をリストに格納して出力します。

 1試行ずつ入れ替え処理を行い、数列を記録します。

# 作図用のオブジェクトを初期化
id_vec          <- 1:N
trace_pivot_vec <- rep(NA, times = N)
trace_df        <- tmp_df

# 試行回数を初期化
iter <- 0

# クイックソート
loop_flag <- TRUE
while(loop_flag) {
  
  # 試行回数をカウント
  iter <- iter + 1
  
  # 前回までのピボット位置を取得
  range_idx_vec <- trace_pivot_vec |> 
    (\(vec) {vec[!is.na(vec)]})() |> # 欠損値を除去
    (\(vec) {c(0, vec, N+1)})() # 範囲計算用の値を追加
  
  # 分割数(前回までのピボット数+1)を取得
  max_cnt <- length(range_idx_vec) - 1
  
  # 前回のピボット値を初期化
  pivot_val_vec <- rep(NA, times = N)
  
  # 分割範囲ごとに処理
  for(pivot_cnt in 1:max_cnt) {
    
    # 分割範囲を取得
    l <- range_idx_vec[pivot_cnt] + 1
    r <- range_idx_vec[pivot_cnt+1] - 1
    
    # 範囲がない(ピボットが隣り合う)なら次の範囲に移行
    if(l > r) next
    
    # ピボットを設定
    res_lt <- pivot_split(value = a, index = id_vec, left = l, right = r)
    a[1:N]      <- res_lt[["value"]] # ピボット前後に整理した数列
    id_vec[1:N] <- res_lt[["index"]] # 元のインデックス
    i     <- res_lt[["new_index"]] # 移動後のピボット位置
    #p_idx <- res_lt[["old_index"]] # 移動前のピボット位置
    
    # ピボットを格納
    pivot_val_vec[i]   <- a[i]
    trace_pivot_vec[i] <- i
  }
  
  # 全要素にピボットが割り当てられたら終了
  if(all(!is.na(trace_pivot_vec))) loop_flag <- FALSE
  
  # 数列とピボットを格納
  tmp_df <- tibble::tibble(
    iteration = iter, 
    id        = id_vec, 
    index     = 1:N, 
    value     = a, 
    pivot_flag       = !is.na(pivot_val_vec), # 今回のピボット
    trace_pivot_flag = !is.na(trace_pivot_vec) # 今回までのピボット
  )
  
  # 数列とピボットを記録
  trace_df <- dplyr::bind_rows(trace_df, tmp_df)
  
  # 途中経過を表示
  #print(paste0("--- iteration:", iter, " ---"))
  #print(a)
}

 各試行でピボットに設定された要素のインデックスを trace_pivot_vec に格納していきます。初期値を全て欠損値 NA にしておき、全ての要素がピボットに設定される(欠損値がなくなる)とループを終了します。
 各試行の始めに、trace_pivot_vec に格納されている値を取り出し(欠損値を取り除き)、頭に 0 と尻に N+1 を追加したベクトルを range_idx_vec とします。
 range_idx_vec の最後を除く値+1(右隣)が分割範囲の左端のインデックス l、最初を除く値-1(左隣)が右端のインデックス r に対応します。よって、range_idx_vec の要素数―1回繰り返し処理を行います。ただし、l, r が隣り合う要素の場合は、処理を行わずに次の範囲に移ります。

 pivot_split() で各範囲のピボットの設定と要素の入れ替えを行います。
 数列 a の初期値のインデックス(1 から N の整数)を id_vec としておき、a の要素の入れ替えに対応するように id_vec の要素も入れ替えます。

 試行ごとに、数列やインデックスをデータフレームに保存します。
 その試行でピボット設定された要素は pivot_flag 列を TRUE、その試行までにピボット設定された要素は trace_pivot_flag 列を TRUE にします。(ここの扱いをもう少し上手くやると作図処理が楽になる気がしますが、力尽きました。)

 初期値のときのコードで、ソートした数列をグラフで確認します。

 昇順に並んでいるのが分かります。

 各試行の数列のグラフを並べて推移を確認します。

 ピボットの位置の描画用のデータフレームを作成します。

# ピボット位置を作成
pivot_idx_df <- trace_df |> 
  dplyr::filter(pivot_flag) # 各試行のピボットデータを抽出
pivot_idx_df
## # A tibble: 50 × 6
##    iteration    id index value pivot_flag trace_pivot_flag
##        <dbl> <int> <int> <dbl> <lgl>      <lgl>           
##  1         1    25    41   1.3 TRUE       TRUE            
##  2         2    50    13  -0.9 TRUE       TRUE            
##  3         2    11    50   5.6 TRUE       TRUE            
##  4         3    18     2  -3.6 TRUE       TRUE            
##  5         3    33    19  -0.2 TRUE       TRUE            
##  6         3    45    45   3   TRUE       TRUE            
##  7         4    15     1  -3.8 TRUE       TRUE            
##  8         4    20     7  -2.3 TRUE       TRUE            
##  9         4    28    16  -0.5 TRUE       TRUE            
## 10         4    14    37   1.1 TRUE       TRUE            
## # ℹ 40 more rows

 各試行においてピボットに設定された要素(pivot_flag 列が TRUE の行)を取り出します。

 分割範囲(入替範囲)の描画用のデータフレームを作成します。

# 分割範囲を作成
range_df <- trace_df |> 
  # 範囲計算用の値を追加
  tibble::add_row(
    tidyr::expand_grid(
      iteration = 0:max(trace_df[["iteration"]]), 
      index     = c(0, N+1)
    ) |> # 全ての組み合わせを作成
      dplyr::mutate(
        pivot_flag       = FALSE, 
        trace_pivot_flag = TRUE
      )
  ) |> 
  dplyr::filter(trace_pivot_flag) |> # ピボット済みデータを抽出
  dplyr::arrange(iteration, index) |> # 値の計算用
  dplyr::group_by(iteration) |> # 値の計算用
  # 作図用の値を計算
  dplyr::mutate(
    iteration = iteration + 1, 
    left   = index + 1, 
    right  = dplyr::lead(index, n = 1) - 1, 
    x      = 0.5 * (left + right), 
    y      = 0.5 * (max(a) + min(c(0, a))), 
    width  = right - left + 1, 
    height = max(a) - min(c(0, a)) + 1
  ) |> 
  dplyr::ungroup() |> 
  dplyr::filter(width >= 1) |> 
  dplyr::select(iteration, left, right, x, y, width, height) |> 
  # ピボットデータを追加
  tibble::add_column(
    pivot_idx_df |> 
      dplyr::select(id, index, value)
  )
range_df
## # A tibble: 50 × 10
##    iteration  left right     x     y width height    id index value
##        <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <int> <int> <dbl>
##  1         1     1    50  25.5   0.9    50   10.4    25    41   1.3
##  2         2     1    40  20.5   0.9    40   10.4    50    13  -0.9
##  3         2    42    50  46     0.9     9   10.4    11    50   5.6
##  4         3     1    12   6.5   0.9    12   10.4    18     2  -3.6
##  5         3    14    40  27     0.9    27   10.4    33    19  -0.2
##  6         3    42    49  45.5   0.9     8   10.4    45    45   3  
##  7         4     1     1   1     0.9     1   10.4    15     1  -3.8
##  8         4     3    12   7.5   0.9    10   10.4    20     7  -2.3
##  9         4    14    18  16     0.9     5   10.4    28    16  -0.5
## 10         4    20    40  30     0.9    21   10.4    14    37   1.1
## # ℹ 40 more rows

 分割位置の計算用に、初期値(0番目の試行)を含めた各試行のインデックスに 0N+1 を追加します。

 各試行までにピボット設定された要素(trace_pivot_flag 列が TRUE の行)を取り出します。抽出されたデータの index 列がピボット設定済みの位置を表します。
 分割範囲の左端の位置lは、0・ピボット位置の右隣なので、index 列+1です。また、右端の位置rは、ピボット位置・N+1の左隣なので、index 列を lead() で1行前にズラした-1です。
 ピボット設定済みの要素が隣り合う場合は、分割範囲が存在しないので、lとr(left 列と right 列)の差が0の行を取り除きます。

 1試行分タイミングを遅らせて(iteration 列に 1 を足して)データを作成します。最後の試行では全ての要素がピボット済み(隣り合う)なので除去されています。
 インデックス列が範囲計算用の値なので削除して、ピボットのインデックスと値の列を追加します。

 ピボットの値の描画用のデータフレームを作成します。

# ピボット値を作成
pivot_val_df <- dplyr::bind_cols(
  pivot_idx_df,   # ピボット位置
  range_df |> # 入替範囲
    dplyr::select(left, right)
) |> 
  # 作図用に値を調整
  dplyr::mutate(
    left = left - 0.5, 
    right = right + 0.5
  )
pivot_val_df
## # A tibble: 50 × 8
##    iteration    id index value pivot_flag trace_pivot_flag  left right
##        <dbl> <int> <int> <dbl> <lgl>      <lgl>            <dbl> <dbl>
##  1         1    25    41   1.3 TRUE       TRUE               0.5  50.5
##  2         2    50    13  -0.9 TRUE       TRUE               0.5  40.5
##  3         2    11    50   5.6 TRUE       TRUE              41.5  50.5
##  4         3    18     2  -3.6 TRUE       TRUE               0.5  12.5
##  5         3    33    19  -0.2 TRUE       TRUE              13.5  40.5
##  6         3    45    45   3   TRUE       TRUE              41.5  49.5
##  7         4    15     1  -3.8 TRUE       TRUE               0.5   1.5
##  8         4    20     7  -2.3 TRUE       TRUE               2.5  12.5
##  9         4    28    16  -0.5 TRUE       TRUE              13.5  18.5
## 10         4    14    37   1.1 TRUE       TRUE              19.5  40.5
## # ℹ 40 more rows

 各試行におけるピボット位置と分割範囲のデータを結合します。分割範囲のインデックス(x軸の値)を枠線の位置(x軸の値)に対応するように調整します。

 試行ごとに数列のグラフを作成します。

# 全試行の数列を作図
ggplot() + 
  geom_bar(data = trace_df, 
           mapping = aes(x = index, y = value, fill = factor(value)), stat = "identity") + # 全ての要素
  geom_tile(data = range_df,
            mapping = aes(x = x, y = y, width = width, height = height, 
                          color = factor(value), fill = factor(value), group = factor(id)),
            alpha = 0.1, linewidth = 0.6, linetype = "dashed", show.legend = FALSE) + # 入替範囲
  geom_vline(data = pivot_idx_df, 
             mapping = aes(xintercept = index, color = factor(value)), 
             linewidth = 0.6, linetype = "dotted", show.legend = FALSE)  + # ピボットの位置
  geom_segment(data = pivot_val_df, 
               mapping = aes(x = left, y = value, xend = right, yend = value), 
               color = "red", linewidth = 0.6, linetype = "twodash", show.legend = FALSE)  + # ピボットまでの上限値
  facet_wrap(iteration ~ ., scales = "free_x", labeller = label_both) + # 試行ごとに分割
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "quick sort", 
       fill = "value", 
       x = "index", y = "value")

 facet_wrap() で、試行(iteration 列)ごとに分割して棒グラフを描画します。
 要素の入れ替えが行われた(分割された)範囲を geom_tile() で描画します。(geom_bar() を使って範囲を描画しようとすると、描画自体はできるのですが、width 引数について警告文が出たので避けました。原因は不明です。)

 初期値を0回目として、各試行における「入れ替え後の数列」をそれぞれ棒グラフ(バー)で表します。
 i回目の試行(i+1枚目のグラフ)において、入れ替えの行われた(分割された)範囲を破線の枠、入れ替え後のピボットの位置を点線の垂直線、ピボットの値を赤色の水平線で示します。枠の色と水平線の高さはピボットの値に対応します。

 各グラフを見ると、枠内の要素(バー)がピボットを境界に「ピボット未満の要素が左側」で「ピボット以上の要素が右側」に配置されているのが分かります。また次のグラフを見ると、ピボットの前後を分割しながら全体がソートされるのが分かります。

 各試行の数列のグラフを切り替えて推移を確認します。

 各試行までのピボットの位置の描画用のデータフレームを作成します。

# ピボットの推移を作成
trace_pivot_df <- dplyr::bind_rows(
  # 各試行のピボットデータ
  trace_df |> 
    dplyr::arrange(iteration, index) |> # フレーム調整用
    dplyr::group_by(id) |> # フレーム調整用
    dplyr::mutate(
      pivot_flag = dplyr::lead(pivot_flag, n = 1, default = FALSE)
    ) |> # (表示フレームの調整用に)ピボット割り当て試行を1つ前に変更
    dplyr::ungroup() |> 
    dplyr::filter(pivot_flag) |> # 各試行のピボットデータを抽出
    dplyr::mutate(
      type = "dashed"
    ), 
  # 過去のピボットデータ
  trace_df |> 
    dplyr::filter(trace_pivot_flag) |> # 過去試行のピボットデータを抽出
    dplyr::mutate(
      type = "dotted"
    )
) |> 
  dplyr::arrange(iteration, index)
trace_pivot_df
## # A tibble: 345 × 7
##    iteration    id index value pivot_flag trace_pivot_flag type  
##        <dbl> <int> <int> <dbl> <lgl>      <lgl>            <chr> 
##  1         0    25    25   1.3 TRUE       FALSE            dashed
##  2         1    50    20  -0.9 TRUE       FALSE            dashed
##  3         1    25    41   1.3 TRUE       TRUE             dotted
##  4         1    11    46   5.6 TRUE       FALSE            dashed
##  5         2    18     6  -3.6 TRUE       FALSE            dashed
##  6         2    50    13  -0.9 TRUE       TRUE             dotted
##  7         2    33    27  -0.2 TRUE       FALSE            dashed
##  8         2    25    41   1.3 FALSE      TRUE             dotted
##  9         2    45    45   3   TRUE       FALSE            dashed
## 10         2    11    50   5.6 TRUE       TRUE             dotted
## # ℹ 335 more rows

 各試行においてピボットに設定された要素(pivot_flag 列が TRUE の行)と、各試行までにピボットに設定された要素(trace_pivot_flag 列が TRUE の行)をそれぞれ取り出して結合します。
 trace_df は、各試行における入れ替え後のデータを格納しています。入れ替え前の(1つ前の試行のフレームに)ピボット位置を表示するために、要素ごとに(id 列でグループ化して)、ピボットの設定フラグを lead() で1つ前の試行(1行前)にズラしてから抽出します。
 各試行のピボットと過去のピボットを描き分けるために、線の種類を指定しておきます。

 分割範囲(入替範囲)の描画用のデータフレームを作成します。

# 分割範囲を作成
range_df <- trace_df |> 
  # 範囲計算用の値を追加
  tibble::add_row(
    tidyr::expand_grid(
      iteration = 0:max(trace_df[["iteration"]]), 
      index     = c(0, N+1)
    ) |> # 全ての組み合わせを作成
      dplyr::mutate(
        pivot_flag       = FALSE, 
        trace_pivot_flag = TRUE
      )
  ) |> 
  dplyr::filter(trace_pivot_flag) |> # ピボット済みデータを抽出
  dplyr::arrange(iteration, index) |> # 値の計算用
  dplyr::group_by(iteration) |> # 値の計算用
  # 作図用の値を計算
  dplyr::mutate(
    left   = index + 1, 
    right  = dplyr::lead(index, n = 1) - 1, 
    x      = 0.5 * (left + right), 
    y      = 0.5 * (max(a) + min(c(0, a))), 
    width  = right - left + 1, 
    height = max(a) - min(c(0, a)) + 1
  ) |> 
  dplyr::ungroup() |> 
  dplyr::filter(width >= 1) |> 
  dplyr::select(iteration, left, right, x, y, width, height) |> 
  # ピボットデータを追加
  tibble::add_column(
    trace_df |> 
      dplyr::filter(pivot_flag) |> # 各試行のピボットデータを抽出
      dplyr::select(id, index, value)
  )
range_df
## # A tibble: 50 × 10
##    iteration  left right     x     y width height    id index value
##        <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <int> <int> <dbl>
##  1         0     1    50  25.5   0.9    50   10.4    25    41   1.3
##  2         1     1    40  20.5   0.9    40   10.4    50    13  -0.9
##  3         1    42    50  46     0.9     9   10.4    11    50   5.6
##  4         2     1    12   6.5   0.9    12   10.4    18     2  -3.6
##  5         2    14    40  27     0.9    27   10.4    33    19  -0.2
##  6         2    42    49  45.5   0.9     8   10.4    45    45   3  
##  7         3     1     1   1     0.9     1   10.4    15     1  -3.8
##  8         3     3    12   7.5   0.9    10   10.4    20     7  -2.3
##  9         3    14    18  16     0.9     5   10.4    28    16  -0.5
## 10         3    20    40  30     0.9    21   10.4    14    37   1.1
## # ℹ 40 more rows

 こちらは、試行回数を調整せずにデータを作成します。

 重複ラベルの描画用のデータフレームを作成します。

# 重複ラベルを作成
dup_label_df <- trace_df |> 
  dplyr::arrange(iteration, id) |> # IDの割当用
  dplyr::group_by(iteration, value) |> # IDの割当用
  dplyr::mutate(
    dup_id = dplyr::row_number(id), # 重複要素にIDを割り当て
    dup_num = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, 
      true = paste0("(", dup_id, ")"), 
      false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup() |> 
  dplyr::arrange(iteration, index)
dup_label_df
## # A tibble: 600 × 9
##    iteration    id index value pivot_flag trace_pivot_flag dup_id dup_num
##        <dbl> <int> <int> <dbl> <lgl>      <lgl>             <int>   <int>
##  1         0     1     1   0   FALSE      FALSE                 1       1
##  2         0     2     2   4.6 FALSE      FALSE                 1       1
##  3         0     3     3  -3   FALSE      FALSE                 1       1
##  4         0     4     4   0.5 FALSE      FALSE                 1       2
##  5         0     5     5  -2.2 FALSE      FALSE                 1       3
##  6         0     6     6   0.2 FALSE      FALSE                 1       3
##  7         0     7     7  -1.5 FALSE      FALSE                 1       1
##  8         0     8     8   1.2 FALSE      FALSE                 1       1
##  9         0     9     9   2.9 FALSE      FALSE                 1       1
## 10         0    10    10  -2.2 FALSE      FALSE                 2       3
## # ℹ 590 more rows
## # ℹ 1 more variable: dup_label <chr>

 重複している値の要素にのみ、それぞれ通し番号のラベルを作成します。

 推移のアニメーションを作成します。

# ソートのアニメーションを作図
graph <- ggplot() + 
  geom_tile(data = range_df,
            mapping = aes(x = x, y = y, width = width, height = height, 
                          color = factor(value), fill = factor(value), group = factor(id)),
            alpha = 0.1, linewidth = 1, linetype = "dashed", show.legend = FALSE) + # 入替範囲
  geom_vline(data = trace_pivot_df, 
             mapping = aes(xintercept = index, 
                           color = factor(value), linetype = type, group = factor(id)), 
             linewidth = 1, show.legend = FALSE)  + # ピボット
  geom_bar(data = trace_df, 
           mapping = aes(x = index, y = value, fill = factor(value), group = factor(id)), 
           stat = "identity") + # 全ての要素
  geom_text(data = trace_df, 
            mapping = aes(x = index, y = 0, label = as.character(value), group = factor(id)), 
            vjust = -0.5, size = 5) + # 要素ラベル
  geom_text(data = dup_label_df, 
            mapping = aes(x = index, y = 0, label = dup_label, group = factor(id)), 
            vjust = 1, size = 4) + # 重複ラベル
  gganimate::transition_states(states = iteration, transition_length = 9, state_length = 1, wrap = FALSE) + # フレーム遷移
  gganimate::ease_aes("cubic-in-out") + # 遷移の緩急
  scale_linetype_identity() + 
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "quick sort", 
       subtitle = "iteration : {next_state}", 
       fill = "value", 
       x = "index", y = "value")

# フレーム数を設定
frame_num <- max(trace_df[["iteration"]]) + 1

# 遷移フレーム数を指定
s <- 20

# gif画像を作成
gganimate::animate(
  plot = graph, 
  nframes = (frame_num + 2)*s, start_pause = s, end_pause = s, fps = 20, 
  width = 1500, height = 1200, 
  renderer = gganimate::gifski_renderer()
)

 transition_states() のフレーム制御の引数 states に試行回数列 iteration を指定して、これまでと同様に作図します。
 animate()plot 引数にグラフオブジェクト、nframes 引数にフレーム数を指定して、gif画像を作成します。
 初期値を含むため試行回数は最大値+1です。transition_states() でアニメーションを作成すると、フレーム間の値を補完して、状態の遷移を描画します。遷移(1試行当たり)のフレーム数を s として、試行回数の s 倍の値を nframes 引数に指定します。

 各試行におけるピボットを破線、過去の試行で設定されたピボット(ソート済みの要素)を点線の水平線で示します。
 分割範囲ごとにピボットが1つ設定され、ピボット未満の要素とピボット以上の要素に分けられ、ピボットの前後で分割する度にソートされるのを確認できます。
 また、重複要素の順序が入れ替わる(安定性・stable性がない)のも分かります。

 ここでの試行回数(iterationの値)は、この可視化方法における試行番号(作図処理上の値)であり他の手法と対応していません。

 この記事では、クイックソートを実装しました。次の記事では、ヒープソートを実装します。

12.6 ヒープソートの実装と可視化

 二分ヒープ(バイナリヒープ・binary heap)を用いるヒープソート(heap sort)を実装します。また、ソートの推移のアニメーションを作成して、アルゴリズムを確認します。
 二分木(バイナリツリー・binary tree)については「二分木を作図したい」、ヒープについては「二分ヒープの実装と可視化」を参照してください。

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(gganimate)

 この記事では、基本的に パッケージ名::関数名() の記法を使っているので、パッケージを読み込む必要はありません。ただし、作図コードについてはパッケージ名を省略しているので、ggplot2 を読み込む必要があります。
 また、ネイティブパイプ演算子 |> を使っています。magrittr パッケージのパイプ演算子 %>% に置き換えても処理できますが、その場合は magrittr を読み込む必要があります。

 「ソートの可視化」において、gganimate パッケージを利用してアニメーション(gif画像)を作成します。ソートの実装自体には使いません。

ソートの実装

 まずは、ヒープソートのアルゴリズムを関数として実装します。
 ヒープ化については「二分ヒープの実装と可視化」を参照してください。

 部分木のヒープ化を実装します。

# 部分木のヒープ化の実装
sub_heapify <- function(vec, i, N) {
  
  # 左側の子のインデックスを計算
  child_idx <- i * 2
  
  # 子がなければ終了
  if(child_idx > N) return(vec)
  
  # 値が大きい方の子のインデックスを設定
  if(child_idx+1 <= N & vec[child_idx+1] > vec[child_idx]) { # (&の左がFALSEだと右が処理されずNAにならないので動く)
    child_idx <- child_idx + 1
  }
  
  # 子が親以下なら終了
  if(vec[child_idx] <= vec[i]) return(vec)
  
  # 子と親を入替
  vec[1:N] <- replace(x = vec, list = c(child_idx, i), values = vec[c(i, child_idx)])
  
  # 子を根としてヒープ化
  vec[1:N] <- sub_heapify(vec, child_idx, N)
  
  # 数列を出力
  return(vec)
}

 数列を vec、入れ替え対象のインデックスを i、数列の要素数を N とします。
 i の初期値として部分木の根のインデックスを受け取ります。
 i 番目の頂点(要素)の2つの子の内、左側の子のインデックスを child_idx として計算します。頂点 \(v_i\) の左側の子 \(v_j\) のインデックスは \(j = 2 i\) で求まります。ただし、\(j \gt N\) のとき、頂点 \(v_i\) に子がない(葉である)ことを意味し、入れ替え不要なので、数列をそのまま出力します。
 \(j + 1 \leq N\) のとき右側の子 \(v_{j+1}\) が存在するので、2つの子の値が大きい方のインデックスを childe_idx とします。
 子の値 vec[child_idx] と親の値 vec[i] を比較して、子が親以下なら入れ替え不要なので数列をそのまま出力(再帰処理を終了)します。子が親より大きければ値を入れ替えて、childe_idx 番目の頂点を根とする部分木を sub_heapify() でヒープ化します。

 sub_heapify() を実行すると内部でさらに(入れ替え対象の頂点が葉になるまで再帰的に) sub_heapify() が実行され、親子の入れ替えが繰り返されて部分木がヒープ化されます。

 アルゴリズムの詳しい解説は、本のcode 12.4(前半)などを参照してください。

 ヒープソートを実装します。

# ヒープソートの実装
heap_sort <- function(vec) {
  
  # 数列の要素数を取得
  N <- length(vec)
  
  # 葉を除くノード番号の最大値を設定
  max_i <- N %/% 2
  
  # 全体をヒープ化(1個目の最大値を探索)
  for(i in max_i:1) {
    
    # i番目の頂点を根とする部分木をヒープ化
    vec[1:N] <- sub_heapify(vec, i, N)
  }
  
  # 全体をソート
  for(i in N:2) { # (後から)
    
    # ヒープの先頭(N-i+1個目の最大値)をヒープの末尾(前からi番目・後からN-i+1番目の要素と入替
    vec[1:N] <- replace(x = vec, list = c(1, i), values = vec[c(i, 1)])
    
    # 前からi-1番目までの要素をヒープ化(N-i+2番目の最大値を探索)
    vec[1:(i-1)] <- sub_heapify(vec[1:(i-1)], 1, i-1)
  }
  
  # 数列を出力
  return(vec)
}

 数列を vec、数列の要素数を N とします。

 まずは、受け取った数列全体をヒープ化します。
 各頂点(要素)を逆順に部分木の根として、sub_heapify() で部分木ごとにヒープ化していきます。ただし、葉(子の無い頂点)はヒープ化不要(ヒープ化済みとみなせる)なので、末尾(N番目)の頂点(要素)の親から処理します。\(N\) 番目の頂点 \(v_N\) の親(葉を除く最後の頂点) \(v_j\) のインデックスは \(j = \lfloor \frac{N}{2} \rfloor\) で求められます。\(\lfloor x \rfloor\) は、床関数(floor())で、\(x\) 以下の最大の整数(小数点以下の切り捨て)を表します。

 続いて、ヒープの性質を利用して数列をソートします。
 ヒープ化した数列は、先頭(1番目)の要素が最大値になります。つまり、ヒープ化は最大値の探索に対応します。「ヒープの最大値(先頭)」を「ヒープの末尾」と入れ替えて(数列全体における後側に移動して)、残りの数列を再度ヒープ化(次の最大値を探索)します。最大値を削除した際のヒープ化は、一度の処理で行えます。
 ヒープ化と最大値の入れ替え(最大値の探索と移動)の操作を交互に繰り返すことで、値が大きい順に探索され、数列の右端(N番目)から順番に並べられます。
 最後の要素(要素が1つになった場合)は、ヒープ化と入れ替えは不要です。

 アルゴリズムの詳しい解説は、本のcode 12.4(前半)などを参照してください。

 実装した関数を試してみます。

 要素数を指定して、数列を作成します。

# 要素数を指定
N <- 50

# 数列を生成
a <- sample(x = 1:(2*N), size = N, replace = TRUE)
a
##  [1] 16  7  9 65 43 35  3 92 74 17 17 97 34  7 65 49  2 78 56 58 43 57 40 76 38
## [26] 21 62 40 28 16 55 48 14 91 75  2 28 84 71 22 50 89  8 48 83 77 77 62 17 40

 乱数生成関数 sample() などで数値ベクトルを生成します。

 ヒープソートにより並べ替えます。

# ソート
heap_sort(a)
##  [1]  2  2  3  7  7  8  9 14 16 16 17 17 17 21 22 28 28 34 35 38 40 40 40 43 43
## [26] 48 48 49 50 55 56 57 58 62 62 65 65 71 74 75 76 77 77 78 83 84 89 91 92 97

 目で確認するのはアレなので、組み込み関数のソート結果と比較して確認します。

# 確認
sum(!(heap_sort(a) == sort(a)))
## [1] 0

 !TRUE, FALSE を反転させて sum() で合計することで、一致しない要素数を得られます。結果が 0 であれば全ての要素が一致したことが分かります。

ソートの可視化

 次は、ヒープソートのアルゴリズムをグラフで確認します。
 親子間のインデックスや座標の計算、バイナリツリーの作成については「二分木を作図したい」を参照してください。

 数列の初期値を作成して、グラフで確認します。

 要素数を指定して、数列を作成します。

# 要素数を指定
N <- 25

# 数列を生成
a <- rnorm(n = N, mean = 0, sd = 1) |> 
  round(digits = 1)
table(a)
## a
## -1.6 -1.5 -1.1   -1 -0.9 -0.8 -0.7 -0.5 -0.3 -0.2 -0.1  0.3  0.6  0.7  0.8  0.9 
##    1    1    2    1    1    1    1    2    1    1    1    1    1    1    3    1 
##  1.2  1.3  1.9  2.2 
##    1    1    1    2

 この例では、(負の値を含めるため)正規乱数生成関数 rnorm() で値を生成して、(重複を含めるため)丸め込み関数 round() で小数第一位にしておきます。

 数列をデータフレームに格納します。

# 要素数を取得
N <- length(a)

# 数列を格納
tmp_df <- tibble::tibble(
  operation = 0,   # 試行回数
  id        = 1:N, # 元のインデックス
  index     = id,  # 各試行のインデックス
  value     = a,   # 要素
  heap_num  = N,   # ヒープの要素数
  target_flag = FALSE, # 入替対象
  sorted_flag = FALSE  # ソート済み要素
)
tmp_df
## # A tibble: 25 × 7
##    operation    id index value heap_num target_flag sorted_flag
##        <dbl> <int> <int> <dbl>    <int> <lgl>       <lgl>      
##  1         0     1     1  -0.5       25 FALSE       FALSE      
##  2         0     2     2   0.3       25 FALSE       FALSE      
##  3         0     3     3  -1         25 FALSE       FALSE      
##  4         0     4     4   1.9       25 FALSE       FALSE      
##  5         0     5     5   2.2       25 FALSE       FALSE      
##  6         0     6     6  -0.5       25 FALSE       FALSE      
##  7         0     7     7   1.2       25 FALSE       FALSE      
##  8         0     8     8   0.8       25 FALSE       FALSE      
##  9         0     9     9  -0.8       25 FALSE       FALSE      
## 10         0    10    10   1.3       25 FALSE       FALSE      
## # ℹ 15 more rows

 作図用に、インデックスなどの値とあわせて数列を格納します。数列のみが与えられている場合は、要素数を N とします。

 数列を棒グラフで確認します。

# 数列を作図
ggplot() + 
  geom_bar(data = tmp_df, 
           mapping = aes(x = index, y = value, fill = factor(value)), stat = "identity") + 
  theme(panel.grid.minor.x = element_blank()) + 
  labs(title = "numerical sequence", 
       subtitle = paste0("operation: ", unique(tmp_df[["operation"]])), 
       fill = "value", 
       x = "index", y = "value")

 この要素(バー)を昇順に並べ替えます。

 バイナツリーでも確認します。

# ノードの座標を作成
d <- 0.6
vertex_df <- tibble::tibble(
  value = a, 
  index = 1:length(a), 
  depth = floor(log2(index)), # 縦方向のノード位置
  col_idx = index - 2^depth + 1, # 深さごとのノード番号
  coord_x = (col_idx * 2 - 1) * 1/(2^depth * 2), # 横方向のノード位置
  label_offset = dplyr::if_else(
    condition = index%%2 == 0, true = 0.5+d, false = 0.5-d
  ) # ラベル位置を左右にズラす
)

# エッジの座標を作成
edge_df <- dplyr::bind_rows(
  # 子ノードの座標
  vertex_df |> 
    dplyr::filter(depth > 0) |> # 根を除去
    dplyr::mutate(
      edge_id   = index, 
      node_type = "childe" # (確認用)
    ), 
  # 親ノードの座標
  vertex_df |> 
    dplyr::filter(depth > 0) |> # 根を除去
    dplyr::mutate(
      edge_id   = index, 
      node_type = "parent", # (確認用)
      depth   = depth - 1, # 縦方向のノード位置
      index   = index %/% 2, 
      col_idx = index - 2^depth + 1, # 深さごとのノード番号
      coord_x = (col_idx * 2 - 1) * 1/(2^depth * 2) # 横方向のノード位置
    )
) |> 
  dplyr::select(!label_offset) |> 
  dplyr::arrange(edge_id, depth)

# ツリーの高さを取得
max_h <- floor(log2(N))

# 二分木を作図
d <- 0.1
ggplot() + 
  geom_path(data = edge_df, 
            mapping = aes(x = coord_x, y = depth, group = edge_id), 
            linewidth = 1) + # 辺
  geom_point(data = vertex_df, 
             mapping = aes(x = coord_x, y = depth), 
             size = 9, shape = "circle filled", fill = "white", stroke = 1) + # 頂点
  geom_text(data = vertex_df, 
            mapping = aes(x = coord_x, y = depth, label = value), 
            size = 3) + # 値ラベル
  geom_text(data = vertex_df, 
            mapping = aes(x = coord_x, y = depth, label = index, hjust = label_offset), 
            size = 3, vjust = -2, color = "green4") + # 位置ラベル
  scale_x_continuous(labels = NULL, name = "") + 
  scale_y_reverse(breaks = 0:max_h, minor_breaks = FALSE) + 
  coord_cartesian(xlim = c(0, 1), ylim = c(max_h+d, -d)) + 
  labs(title = "binary tree", 
       subtitle = paste0("operation: ", unique(tmp_df[["operation"]])), 
       y = "depth")

 ヒープ条件(それぞれの分岐で親が子より大きい)を満たしていません。

 ヒープ化(数列の入れ替え)時に、値と共にインデックスも入れ替えて出力する関数を定義しておきます。

# 部分木のヒープ化の実装
trace_sub_heapify <- function(value, index, i, N) {
  
  # 値・インデックスを取得
  val_vec <- value
  idx_vec <- index
  
  # 左側の子のインデックスを計算
  child_idx <- i * 2
  
  # 子がなければ終了
  if(child_idx > N) {
    lt <- list(value = val_vec, index = idx_vec)
    return(lt)
  }
  
  # 値が大きい方の子のインデックスを設定
  if(child_idx+1 <= N & val_vec[child_idx+1] > val_vec[child_idx]) { # (&の左がFALSEだと右が処理されずNAにならないので動く)
    child_idx <- child_idx + 1
  }
  
  # 子が親以下なら終了
  if(val_vec[child_idx] <= val_vec[i]) {
    lt <- list(value = val_vec, index = idx_vec)
    return(lt)
  }
  
  # 子と親を入替
  val_vec[1:N] <- replace(x = val_vec, list = c(child_idx, i), values = val_vec[c(i, child_idx)])
  idx_vec[1:N] <- replace(x = idx_vec, list = c(child_idx, i), values = idx_vec[c(i, child_idx)])
  
  # 子を根としてヒープ化
  res_lt <- trace_sub_heapify(val_vec, idx_vec, child_idx, N)
  val_vec[1:N] <- res_lt[["value"]]
  idx_vec[1:N] <- res_lt[["index"]]
  
  # 数列を出力
  lt <- list(value = val_vec, index = idx_vec)
  return(lt)
}
# ヒープ化の実装
trace_heapify <- function(value, index) {
  
  # 値・インデックスを取得
  val_vec <- value
  idx_vec <- index
  N       <- length(val_vec)
  
  # 葉を除くノード番号の最大値を設定
  max_i <- N %/% 2
  
  # 全体をヒープ化
  for(i in max_i:1) {
    
    # i番目の頂点を根とする部分木をヒープ化
    res_lt <- trace_sub_heapify(val_vec, idx_vec, i, N)
    val_vec[1:N] <- res_lt[["value"]]
    idx_vec[1:N] <- res_lt[["index"]]
  }
  
  # 数列を出力
  lt <- list(value = val_vec, index = idx_vec)
  return(lt)
}

 数列の値とインデックスを受け取り、sub_heap()heap_sort() の前半と同様に入れ替えを行い、作図用の値を出力します。
 作図用の値として、「数列の値・インデックス」のベクトルをリストに格納して出力します。

 アルゴリズムに従いソートを行って、作図用のデータを作成します。

 数列全体をヒープ化します。

# 作図用のオブジェクトを初期化
id_vec   <- 1:N
trace_df <- tmp_df

# 数列をヒープ化
res_lt <- trace_heapify(a, id_vec)
a[1:N]      <- res_lt[["value"]]
id_vec[1:N] <- res_lt[["index"]]

# 数列を格納
tmp_df <- tibble::tibble(
  operation = 1, 
  id        = id_vec, 
  index     = 1:N, 
  value     = a, 
  heap_num  = N, 
  target_flag = FALSE, 
  sorted_flag = FALSE
)
tmp_df
## # A tibble: 25 × 7
##    operation    id index value heap_num target_flag sorted_flag
##        <dbl> <int> <int> <dbl>    <int> <lgl>       <lgl>      
##  1         1     5     1   2.2       25 FALSE       FALSE      
##  2         1    20     2   2.2       25 FALSE       FALSE      
##  3         1     7     3   1.2       25 FALSE       FALSE      
##  4         1     4     4   1.9       25 FALSE       FALSE      
##  5         1    10     5   1.3       25 FALSE       FALSE      
##  6         1    13     6   0.8       25 FALSE       FALSE      
##  7         1    15     7   0.9       25 FALSE       FALSE      
##  8         1     8     8   0.8       25 FALSE       FALSE      
##  9         1    19     9   0.8       25 FALSE       FALSE      
## 10         1    21    10   0.7       25 FALSE       FALSE      
## # ℹ 15 more rows

 数列のヒープ化に対応するようにインデックスも入れ替えて、データフレームに格納します。

 先ほどのコードで、数列を棒グラフとツリーで確認します。

 数列全体がヒープ化された状態です。

 1試行ずつ入れ替え処理を行い、数列を記録します。

# 数列を記録
trace_df <- dplyr::bind_rows(trace_df, tmp_df)

# 全体をソート
iter <- 1
for(i in N:2) {
  
  ## 最大値の入替操作
  
  # N-i+1個目の最大値をヒープの末尾と入替
  a[1:N]      <- replace(x = a, list = c(1, i), values = a[c(i, 1)])
  id_vec[1:N] <- replace(x = id_vec, list = c(1, i), values = id_vec[c(i, 1)])
  
  # 数列を格納
  tmp_df <- tibble::tibble(
    operation = iter + 1, 
    id        = id_vec, 
    index     = 1:N, 
    value     = a, 
    heap_num  = i, # (作図用に入れ替え後もi番目をヒープとしておく)
    target_flag = index %in% c(1, i), # 入替要素:(最大値と末尾)
    sorted_flag = index >= i
  )
  
  # 数列を記録
  trace_df <- dplyr::bind_rows(trace_df, tmp_df)
  
  # 途中経過を表示
  #print(paste0("--- operation: ", iter+1, " ---"))
  #print(a)
  
  ## ヒープ化(最大値の探索)操作
  
  # 要素が1つならヒープ化不要
  if(i == 2) break
  
  # i-1番目までの要素をヒープ化(N-i+2個目の最大値を探索)
  res_lt <- trace_sub_heapify(a[1:(i-1)], id_vec[1:(i-1)], i = 1, N = i-1)
  a[1:(i-1)]      <- res_lt[["value"]]
  id_vec[1:(i-1)] <- res_lt[["index"]]
  
  # 数列を格納
  tmp_df <- tibble::tibble(
    operation = iter + 2, 
    id        = id_vec, 
    index     = 1:N, 
    value     = a, 
    heap_num  = i - 1, 
    target_flag = FALSE, 
    sorted_flag = index >= i
  )
  
  # 数列を記録
  trace_df <- dplyr::bind_rows(trace_df, tmp_df)
  
  # 試行回数をカウント
  iter <- iter + 2
}

 「ソートの実装」(heap_sort() の後半)のときと同様にして、要素を入れ替えていきます。
 数列 a の初期値のインデックス(1 から N の整数)を id_vec としておき、a の要素の入れ替えに対応するように id_vec の要素も入れ替えます。

 試行ごとに、数列やインデックスをデータフレームに保存します。
 各試行における、ヒープ部分の先頭(最大値)と末尾の要素は target_flag 列を TRUE、ソートされた(ヒープ部分でない)要素は sorted_flag 列を TRUE にします。

 初期値のときのコードで、最終結果の数列を棒グラフとツリーで確認します。

 数列全体が昇順にソートされた状態です。ヒープ条件は満たしていません。

 各試行の数列の棒グラフを切り替えて推移を確認します。

 重複ラベルの描画用のデータフレームを作成します。

# 重複ラベルを作成
dup_label_df <- trace_df |> 
  dplyr::arrange(operation, id) |> # IDの割当用
  dplyr::group_by(operation, value) |> # IDの割当用
  dplyr::mutate(
    dup_id    = dplyr::row_number(id), # 重複IDを割当
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup() |> 
  dplyr::arrange(operation, index)
dup_label_df |> 
  dplyr::select(!heap_num) # (資料作成用に間引き)
## # A tibble: 1,225 × 9
##    operation    id index value target_flag sorted_flag dup_id dup_num dup_label
##        <dbl> <int> <int> <dbl> <lgl>       <lgl>        <int>   <int> <chr>    
##  1         0     1     1  -0.5 FALSE       FALSE            1       2 "(1)"    
##  2         0     2     2   0.3 FALSE       FALSE            1       1 ""       
##  3         0     3     3  -1   FALSE       FALSE            1       1 ""       
##  4         0     4     4   1.9 FALSE       FALSE            1       1 ""       
##  5         0     5     5   2.2 FALSE       FALSE            1       2 "(1)"    
##  6         0     6     6  -0.5 FALSE       FALSE            2       2 "(2)"    
##  7         0     7     7   1.2 FALSE       FALSE            1       1 ""       
##  8         0     8     8   0.8 FALSE       FALSE            1       3 "(1)"    
##  9         0     9     9  -0.8 FALSE       FALSE            1       1 ""       
## 10         0    10    10   1.3 FALSE       FALSE            1       1 ""       
## # ℹ 1,215 more rows

 各試行の数列データ trace_df を用いて、重複している値の要素にのみそれぞれ通し番号のラベルを作成します。

 最大値と入替対象の描画用のデータフレームを作成します。

# 入替対象の座標を作成
range_swap_df <- trace_df |> 
  dplyr::filter(operation%%2 == 1) |> # 入替前のデータを抽出
  dplyr::group_by(operation) |> # 入替対象の抽出用
  dplyr::filter(index %in% c(1, N-0.5*(operation-1))) |> # 入替対象の抽出
  dplyr::ungroup() |> 
  dplyr::mutate(
    bar_label = paste0("iter: ", operation+1, ", idx: ", index)
  )
range_swap_df
## # A tibble: 48 × 8
##    operation    id index value heap_num target_flag sorted_flag bar_label       
##        <dbl> <int> <int> <dbl>    <dbl> <lgl>       <lgl>       <chr>           
##  1         1     5     1   2.2       25 FALSE       FALSE       iter: 2, idx: 1 
##  2         1    25    25  -1.5       25 FALSE       FALSE       iter: 2, idx: 25
##  3         3    20     1   2.2       24 FALSE       FALSE       iter: 4, idx: 1 
##  4         3    24    24  -1.1       24 FALSE       FALSE       iter: 4, idx: 24
##  5         5     4     1   1.9       23 FALSE       FALSE       iter: 6, idx: 1 
##  6         5    23    23  -1.6       23 FALSE       FALSE       iter: 6, idx: 23
##  7         7    10     1   1.3       22 FALSE       FALSE       iter: 8, idx: 1 
##  8         7    22    22  -0.9       22 FALSE       FALSE       iter: 8, idx: 22
##  9         9     7     1   1.2       21 FALSE       FALSE       iter: 10, idx: 1
## 10         9     1    21  -0.5       21 FALSE       FALSE       iter: 10, idx: …
## # ℹ 38 more rows

 各試行における入れ替え前(ヒープ化後)の数列全体(operation 列が奇数)の行を取り出して、ヒープの先頭と末尾(index 列が \(N - \frac{i - 1}{2}\) )の行を取り出します。(target_flag 列が TRUE なのは入れ替え後のデータです。)

 最大値の描画用のデータフレームを作成します。

# ヒープの最大値を作成
upper_df <- trace_df |> 
  dplyr::filter(target_flag) |> # 入替時のデータを抽出
  dplyr::group_by(operation) |> # 最大値の抽出用
  dplyr::filter(index != 1) |> # 最大値を抽出
  dplyr::ungroup() |> 
  dplyr::mutate(
    operation = operation - 1 # 表示フレームを調整
  )
upper_df
## # A tibble: 24 × 7
##    operation    id index value heap_num target_flag sorted_flag
##        <dbl> <int> <int> <dbl>    <dbl> <lgl>       <lgl>      
##  1         1     5    25   2.2       25 TRUE        TRUE       
##  2         3    20    24   2.2       24 TRUE        TRUE       
##  3         5     4    23   1.9       23 TRUE        TRUE       
##  4         7    10    22   1.3       22 TRUE        TRUE       
##  5         9     7    21   1.2       21 TRUE        TRUE       
##  6        11    15    20   0.9       20 TRUE        TRUE       
##  7        13     8    19   0.8       19 TRUE        TRUE       
##  8        15    19    18   0.8       18 TRUE        TRUE       
##  9        17    13    17   0.8       17 TRUE        TRUE       
## 10        19    21    16   0.7       16 TRUE        TRUE       
## # ℹ 14 more rows

 入れ替え後(target_flag 列が TRUE)のデータ(行)を取り出して、最大値(index 列が 1 でない)の行を取り出します。

 ヒープ化要素の描画用のデータフレームを作成します。

# ヒープ範囲の座標を作成
d <- 0.2
range_heap_df <- tibble::tibble(
  operation = 0:(2*(N-1)), 
  left  = 1, 
  right = c(NA, N, rep((N-1):2, each = 2), 0), 
  x     = 0.5 * (left + right), 
  y     = 0.5 * (max(a) + min(c(0, a))), 
  w     = right - left + 1, 
  h     = max(a) - min(c(0, a)) + d
)
range_heap_df
## # A tibble: 49 × 7
##    operation  left right     x     y     w     h
##        <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1         0     1    NA  NA     0.3    NA     4
##  2         1     1    25  13     0.3    25     4
##  3         2     1    24  12.5   0.3    24     4
##  4         3     1    24  12.5   0.3    24     4
##  5         4     1    23  12     0.3    23     4
##  6         5     1    23  12     0.3    23     4
##  7         6     1    22  11.5   0.3    22     4
##  8         7     1    22  11.5   0.3    22     4
##  9         8     1    21  11     0.3    21     4
## 10         9     1    21  11     0.3    21     4
## # ℹ 39 more rows

 ヒープの要素数は、初回時(1回目・operation 列が 1)はN個であり、入替と探索の2回ごとに1個ずつ減っていきます。ただし、初期値(0回目)は0個で、最終試行(2N-1回目)では2個減り0個になります。
 各試行におけるヒープ要素(バー)を囲うように座標を作成します。

 ソート済み要素の描画用のデータフレームを作成します。

# ソート済み範囲の座標を作成
range_sort_df <- tibble::tibble(
  operation = 0:(2*(N-1)), 
  left  = c(NA, N+1, rep(N:3, each = 2), 1), 
  right = N, 
  x     = 0.5 * (left + right), 
  y     = 0.5 * (max(a) + min(c(0, a))), 
  w     = right - left + 1, 
  h     = max(a) - min(c(0, a)) + d
)
range_sort_df
## # A tibble: 49 × 7
##    operation  left right     x     y     w     h
##        <int> <dbl> <int> <dbl> <dbl> <dbl> <dbl>
##  1         0    NA    25  NA     0.3    NA     4
##  2         1    26    25  25.5   0.3     0     4
##  3         2    25    25  25     0.3     1     4
##  4         3    25    25  25     0.3     1     4
##  5         4    24    25  24.5   0.3     2     4
##  6         5    24    25  24.5   0.3     2     4
##  7         6    23    25  24     0.3     3     4
##  8         7    23    25  24     0.3     3     4
##  9         8    22    25  23.5   0.3     4     4
## 10         9    22    25  23.5   0.3     4     4
## # ℹ 39 more rows

 ヒープの要素数とは逆に右端から数えて、ソート済み要素数は、初回時(1回目・operation 列が 1)は0個であり、入替と探索の2回ごとに1個ずつ増えていきます。ただし、初期値(0回目)は0個で、最終試行(2N-1回目)では2個増えN個になります。
 各試行におけるソート済み要素(バー)を囲うように座標を作成します。

 推移のアニメーションを作成します。

# ソートのアニメーションを作図
graph <- ggplot() + 
  geom_tile(data = range_heap_df,
            mapping = aes(x = x, y = y, width = w, height = h, color = "heap"), 
            fill = "blue", alpha = 0.1, linewidth = 1, linetype = "dashed", na.rm = TRUE) + # ヒープ範囲
  geom_tile(data = range_sort_df,
            mapping = aes(x = x, y = y, width = w, height = h, color = "sort"), 
            fill = "orange", alpha = 0.1, linewidth = 1, linetype = "dashed", na.rm = TRUE) + # ソート済み範囲
  geom_bar(data = range_swap_df,
           mapping = aes(x = index, y = value, color = "max", group = bar_label), 
           stat = "identity", 
           fill = "red", alpha = 0.1, linewidth = 1, linetype = "dashed") + # 入替対象
  geom_bar(data = trace_df,
           mapping = aes(x = index, y = value, fill = factor(value), group = id),
           stat = "identity", show.legend = FALSE) + # 全ての要素
  geom_text(data = trace_df,
            mapping = aes(x = index, y = 0, label = as.character(value), group = id),
            vjust = -0.5, size = 4) + # 要素ラベル
  geom_text(data = dup_label_df,
            mapping = aes(x = index, y = 0, label = dup_label, group = id),
            vjust = 1, size = 3) + # 重複ラベル
  geom_hline(data = upper_df, 
             mapping = aes(yintercept = value), 
             color = "red", linewidth = 1, linetype = "dotted") + # ヒープの最大値
  gganimate::transition_states(states = operation, transition_length = 9, state_length = 1, wrap = FALSE) + # フレーム遷移
  gganimate::ease_aes("cubic-in-out") + # 遷移の緩急
  scale_color_manual(breaks = c("max", "heap", "sort"), 
                     values = c("red", "blue", "orange"), 
                     labels = c("swap", "heap", "sorted"), 
                     name = "bar") + # (凡例表示用)
  guides(color = guide_legend(override.aes = list(fill = c("red", "blue", "orange")))) + # 凡例の体裁
  theme(panel.grid.minor.x = element_blank()) + 
  labs(title = "heap sort", 
       subtitle = "operation: {next_state}", 
       x = "index", y = "value")

# フレーム数を設定
frame_num <- max(trace_df[["operation"]]) + 1
frame_num <- 2 * N - 1

# 遷移フレーム数を指定
s <- 20

# gif画像を作成
gganimate::animate(
  plot = graph, 
  nframes = (frame_num + 2)*s, start_pause = s, end_pause = s, fps = 20, 
  width = 900, height = 600, 
  renderer = gganimate::gifski_renderer()
)

 transition_states() のフレーム制御の引数 states に試行回数列 operation を指定して、これまでと同様に作図します。
 animate()plot 引数にグラフオブジェクト、nframes 引数にフレーム数を指定して、gif画像を作成します。
 初期値を含み・要素ごとに最大値の探索と入替の繰り返し・ただし1要素の場合は処理が不要なため、試行回数は \(1 + 2 (N - 1) = 2 N - 1\) です。transition_states() でアニメーションを作成すると、フレーム間の値を補完して、状態の遷移を描画します。遷移のフレームを s として、試行回数の s 倍の値を nframes 引数に指定します。

 ヒープとして扱う要素を青色、ヒープの先頭(最大値)と末尾の要素を赤色、ソート済みの要素をオレンジ色の塗りつぶし枠で示します。また、各試行における最大値を赤色の水平線で示します。
 最大値の探索(ヒープ化)と入替(削除)を繰り返すことで、降順に(右端から)ソートされるのを確認できます。青色の枠内のバー(要素)は赤線よりも小さく、オレンジ色の枠内のバーよりも大きいのが分かります。
 また、重複要素の順序が入れ替わる(安定性・stable性がない)も確認できます。

 ここでの試行回数(operationの値)は、この可視化方法における試行番号(作図処理上の値)であり他の手法と対応していません。

 続いて、各試行の数列のバイナツリーを切り替えて推移を確認します。

 試行ごとの頂点の描画用のデータフレームを作成します。

# 頂点の座標を作成
trace_vertex_df <- trace_df |> 
  dplyr::mutate(
    depth   = floor(log2(index)), # 縦方向のノード位置
    col_idx = index - 2^depth + 1, 
    coord_x = (col_idx * 2 - 1) * 1/(2^depth * 2), # 横方向のノード位置
  ) |> 
  dplyr::arrange(operation, index)
trace_vertex_df |> 
  dplyr::select(!heap_num) # (資料作成用に間引き)
## # A tibble: 1,225 × 9
##    operation    id index value target_flag sorted_flag depth col_idx coord_x
##        <dbl> <int> <int> <dbl> <lgl>       <lgl>       <dbl>   <dbl>   <dbl>
##  1         0     1     1  -0.5 FALSE       FALSE           0       1  0.5   
##  2         0     2     2   0.3 FALSE       FALSE           1       1  0.25  
##  3         0     3     3  -1   FALSE       FALSE           1       2  0.75  
##  4         0     4     4   1.9 FALSE       FALSE           2       1  0.125 
##  5         0     5     5   2.2 FALSE       FALSE           2       2  0.375 
##  6         0     6     6  -0.5 FALSE       FALSE           2       3  0.625 
##  7         0     7     7   1.2 FALSE       FALSE           2       4  0.875 
##  8         0     8     8   0.8 FALSE       FALSE           3       1  0.0625
##  9         0     9     9  -0.8 FALSE       FALSE           3       2  0.188 
## 10         0    10    10   1.3 FALSE       FALSE           3       3  0.312 
## # ℹ 1,215 more rows

 各試行の数列データ trace_df を用いて、試行ごとに(step 列でグループ化して)、「ヒープの可視化」のときと同様にして、頂点の座標を計算します。

 辺の描画用のデータフレームを作成します。

# 辺の座標を作成
edge_df <- dplyr::bind_rows(
  # 子ノードの座標
  trace_vertex_df |> 
    dplyr::filter(operation == 0, depth > 0) |> # 1試行分のデータを抽出・根を除去
    dplyr::mutate(
      edge_id = index
    ), 
  # 親ノードの座標
  trace_vertex_df |> 
    dplyr::filter(operation == 0, depth > 0) |> # 1試行分のデータを抽出・根を除去
    dplyr::mutate(
      edge_id = index, # 子インデックス
      depth   = depth - 1, 
      index   = index %/% 2, # 親インデックス
      col_idx = index - 2^depth + 1, 
      coord_x = (col_idx * 2 - 1) * 1/(2^depth * 2)
    )
) |> 
  dplyr::select(!operation) |> # フレーム遷移の影響から外す
  dplyr::arrange(edge_id, depth)
edge_df |> 
  dplyr::select(!col_idx) # (資料作成用に間引き)
## # A tibble: 48 × 9
##       id index value heap_num target_flag sorted_flag depth coord_x edge_id
##    <int> <dbl> <dbl>    <dbl> <lgl>       <lgl>       <dbl>   <dbl>   <int>
##  1     2     1   0.3       25 FALSE       FALSE           0   0.5         2
##  2     2     2   0.3       25 FALSE       FALSE           1   0.25        2
##  3     3     1  -1         25 FALSE       FALSE           0   0.5         3
##  4     3     3  -1         25 FALSE       FALSE           1   0.75        3
##  5     4     2   1.9       25 FALSE       FALSE           1   0.25        4
##  6     4     4   1.9       25 FALSE       FALSE           2   0.125       4
##  7     5     2   2.2       25 FALSE       FALSE           1   0.25        5
##  8     5     5   2.2       25 FALSE       FALSE           2   0.375       5
##  9     6     3  -0.5       25 FALSE       FALSE           1   0.75        6
## 10     6     6  -0.5       25 FALSE       FALSE           2   0.625       6
## # ℹ 38 more rows

 1試行分(この例だと初期値)の頂点の座標データを取り出して、「ヒープの可視化」のときと同様にして、辺(始点・終点)の座標を計算します。

 インデックスラベルの描画用のデータフレームを作成します。

# インデックスラベルの座標を作成
d <- 0.6
index_df <- trace_vertex_df |> 
  dplyr::filter(operation == 0) |> # 1試行分のデータを抽出
  dplyr::mutate(
    offset = dplyr::if_else(
      condition = index%%2 == 0, true = 0.5+d, false = 0.5-d
    ) # ラベル位置を左右にズラす
  ) |> 
  dplyr::select(!operation) # フレーム遷移の影響から外す
index_df |> 
  dplyr::select(!col_idx) # (資料作成用に間引き)
## # A tibble: 25 × 9
##       id index value heap_num target_flag sorted_flag depth coord_x offset
##    <int> <int> <dbl>    <dbl> <lgl>       <lgl>       <dbl>   <dbl>  <dbl>
##  1     1     1  -0.5       25 FALSE       FALSE           0  0.5      -0.1
##  2     2     2   0.3       25 FALSE       FALSE           1  0.25      1.1
##  3     3     3  -1         25 FALSE       FALSE           1  0.75     -0.1
##  4     4     4   1.9       25 FALSE       FALSE           2  0.125     1.1
##  5     5     5   2.2       25 FALSE       FALSE           2  0.375    -0.1
##  6     6     6  -0.5       25 FALSE       FALSE           2  0.625     1.1
##  7     7     7   1.2       25 FALSE       FALSE           2  0.875    -0.1
##  8     8     8   0.8       25 FALSE       FALSE           3  0.0625    1.1
##  9     9     9  -0.8       25 FALSE       FALSE           3  0.188    -0.1
## 10    10    10   1.3       25 FALSE       FALSE           3  0.312     1.1
## # ℹ 15 more rows

 1試行分の頂点の座標データを取り出して、「ヒープの可視化」のときと同様にして、ラベル位置の調整値を追加します。

 重複ラベルの描画用のデータフレームを作成します。

# 重複ラベルを作成
dup_label_df <- trace_vertex_df |> 
  dplyr::arrange(operation, id) |> # IDの割当用
  dplyr::group_by(operation, value) |> # IDの割当用
  dplyr::mutate(
    dup_id    = dplyr::row_number(id), # 重複IDを割当
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup() |> 
  dplyr::arrange(operation, index)
dup_label_df |> 
  dplyr::select(!c(target_flag, sorted_flag, col_idx, coord_x)) # (資料作成用に間引き)
## # A tibble: 1,225 × 9
##    operation    id index value heap_num depth dup_id dup_num dup_label
##        <dbl> <int> <int> <dbl>    <dbl> <dbl>  <int>   <int> <chr>    
##  1         0     1     1  -0.5       25     0      1       2 "(1)"    
##  2         0     2     2   0.3       25     1      1       1 ""       
##  3         0     3     3  -1         25     1      1       1 ""       
##  4         0     4     4   1.9       25     2      1       1 ""       
##  5         0     5     5   2.2       25     2      1       2 "(1)"    
##  6         0     6     6  -0.5       25     2      2       2 "(2)"    
##  7         0     7     7   1.2       25     2      1       1 ""       
##  8         0     8     8   0.8       25     3      1       3 "(1)"    
##  9         0     9     9  -0.8       25     3      1       1 ""       
## 10         0    10    10   1.3       25     3      1       1 ""       
## # ℹ 1,215 more rows

 「棒グラフによる可視化」のときと同様にして、重複ラベルを作成します。

 入れ替え対象・ヒープ化済み・ソート済みの頂点の描画用のデータフレームを作成します。

# 入替対象の頂点の座標を作成
target_vertex_df <- dplyr::bind_rows(
  # 最大値との入替要素の座標
  trace_vertex_df |> 
    dplyr::filter(target_flag) |> # 入替要素を抽出
    dplyr::mutate(
      vertex_type = "max"
    ), 
  # (ヒープの最大値・末尾を除く)ヒープ化要素の座標
  trace_vertex_df |> 
    dplyr::filter(operation > 1, !sorted_flag, !target_flag) |> # ヒープ要素を抽出
    dplyr::mutate(
      vertex_type = "heap"
    ), 
  # ソート済み要素の座標
  trace_vertex_df |> 
    dplyr::filter(sorted_flag, !target_flag) |> # ソート済み要素を抽出
    dplyr::mutate(
      vertex_type = "sort"
    ), 
  # フレーム調整用に最終値を複製
  trace_vertex_df |> 
    dplyr::filter(operation == max(operation)) |> # 最終結果を抽出
    dplyr::mutate(
      vertex_type = "sort", 
      heap_num    = 0, 
      target_flag = FALSE, 
      sorted_flag = TRUE, 
      operation   = operation + 1 # フレーム調整用
    ), 
  tibble::tibble(
    operation = 1
  ) # (フレーム順の謎バグ回避用)
) |> 
  dplyr::mutate(
    vertex_label = paste0("iter: ", operation, ", idx: ", index), 
    operation = operation - 1 # 表示フレームを調整
  ) |> 
  dplyr::arrange(operation, index)
target_vertex_df |> 
  dplyr::select(!c(value, heap_num, col_idx, coord_x)) # (資料作成用に間引き)
## # A tibble: 1,201 × 8
##    operation    id index target_flag sorted_flag depth vertex_type vertex_label 
##        <dbl> <int> <int> <lgl>       <lgl>       <dbl> <chr>       <chr>        
##  1         0    NA    NA NA          NA             NA <NA>        iter: 1, idx…
##  2         1    25     1 TRUE        FALSE           0 max         iter: 2, idx…
##  3         1    20     2 FALSE       FALSE           1 heap        iter: 2, idx…
##  4         1     7     3 FALSE       FALSE           1 heap        iter: 2, idx…
##  5         1     4     4 FALSE       FALSE           2 heap        iter: 2, idx…
##  6         1    10     5 FALSE       FALSE           2 heap        iter: 2, idx…
##  7         1    13     6 FALSE       FALSE           2 heap        iter: 2, idx…
##  8         1    15     7 FALSE       FALSE           2 heap        iter: 2, idx…
##  9         1     8     8 FALSE       FALSE           3 heap        iter: 2, idx…
## 10         1    19     9 FALSE       FALSE           3 heap        iter: 2, idx…
## # ℹ 1,191 more rows

 各試行の頂点の座標データ trace_vertex_df から、最大値との入れ替え要素・ヒープ要素・ソート済み要素をそれぞれ取り出して結合します。
 入れ替え要素は、target_flag 列が TRUE の行を取り出します。
 ヒープ要素は、sorted_flag 列が FALSE の行を取り出します。ただし、初期値と初回時のデータ、各試行のヒープの最大値は除きます。
 ソート済み要素は、sorted_flag 列が TRUE の行を取り出します。
 1試行早く表示するために、step 列を-1し、最終値(全ての要素がソート済み)のデータを複製します。

 入れ替え対象の頂点を結ぶ辺の描画用のデータフレームを作成します。

# 最大値との入替対象の辺の座標を計算
swap_edge_df <- target_vertex_df |> 
  dplyr::filter(target_flag) |> # 入替要素を抽出
  dplyr::arrange(operation, index)
swap_edge_df |> 
  dplyr::select(!c(value, heap_num, col_idx, coord_x)) # (資料作成用に間引き)
## # A tibble: 48 × 8
##    operation    id index target_flag sorted_flag depth vertex_type vertex_label 
##        <dbl> <int> <int> <lgl>       <lgl>       <dbl> <chr>       <chr>        
##  1         1    25     1 TRUE        FALSE           0 max         iter: 2, idx…
##  2         1     5    25 TRUE        TRUE            4 max         iter: 2, idx…
##  3         3    24     1 TRUE        FALSE           0 max         iter: 4, idx…
##  4         3    20    24 TRUE        TRUE            4 max         iter: 4, idx…
##  5         5    23     1 TRUE        FALSE           0 max         iter: 6, idx…
##  6         5     4    23 TRUE        TRUE            4 max         iter: 6, idx…
##  7         7    22     1 TRUE        FALSE           0 max         iter: 8, idx…
##  8         7    10    22 TRUE        TRUE            4 max         iter: 8, idx…
##  9         9     1     1 TRUE        FALSE           0 max         iter: 10, id…
## 10         9     7    21 TRUE        TRUE            4 max         iter: 10, id…
## # ℹ 38 more rows

 各試行の頂点ごとの属性・座標データ target_vertex_df から、各試行の入れ替え要素(target_flag 列が TRUE)の行を取り出します。(これは試行回数列が欠損していても問題ない理由が分からない。)

 ヒープ化済みの頂点を結ぶ辺の描画用のデータフレームを作成します。

# ヒープ対象の辺の座標を計算
heap_edge_df <- dplyr::bind_rows(
  # 親の座標
  target_vertex_df |> 
    dplyr::filter(index <= heap_num%/%2) |> # 葉とソート済み要素を除去
    tidyr::uncount(weights = 2, .id = "childe_id") |> # 子の数に複製
    dplyr::mutate(
      vertex_type = "heap", 
      childe_id = childe_id - 1, # 行番号を子IDに変換
      edge_id   = index*2+childe_id # 子インデックス
    ) |> 
  dplyr::filter(edge_id <= heap_num), # 子がなければ除去
  # 子の座標
  target_vertex_df |> 
    dplyr::filter(index > 1, index <= heap_num) |> # 根とソート済み要素を除去
    dplyr::mutate(
      vertex_type = "heap", 
      edge_id = index
    ), 
  tibble::tibble(
    operation = 0
  ) # (フレーム順の謎バグ回避用)
) |> 
  dplyr::mutate(
    edge_label = paste0("iter: ", operation+1, ", idx: ", edge_id)
  ) |> 
  dplyr::arrange(operation, edge_id, index)
heap_edge_df |> 
  dplyr::select(!c(value, heap_num, depth, col_idx, coord_x, vertex_type, vertex_label)) # (資料作成用に間引き)
## # A tibble: 1,153 × 8
##    operation    id index target_flag sorted_flag childe_id edge_id edge_label   
##        <dbl> <int> <int> <lgl>       <lgl>           <dbl>   <dbl> <chr>        
##  1         0    NA    NA NA          NA                 NA      NA iter: 1, idx…
##  2         1    25     1 TRUE        FALSE               0       2 iter: 2, idx…
##  3         1    20     2 FALSE       FALSE              NA       2 iter: 2, idx…
##  4         1    25     1 TRUE        FALSE               1       3 iter: 2, idx…
##  5         1     7     3 FALSE       FALSE              NA       3 iter: 2, idx…
##  6         1    20     2 FALSE       FALSE               0       4 iter: 2, idx…
##  7         1     4     4 FALSE       FALSE              NA       4 iter: 2, idx…
##  8         1    20     2 FALSE       FALSE               1       5 iter: 2, idx…
##  9         1    10     5 FALSE       FALSE              NA       5 iter: 2, idx…
## 10         1     7     3 FALSE       FALSE               0       6 iter: 2, idx…
## # ℹ 1,143 more rows

 ヒープ要素の親と子のデータをそれぞれ取り出して結合します。sorted_flag 列が FALSE のデータは入れ替え後の最大値を含みません。
 親の頂点として「1番目から葉を除いた(末尾の親までの)要素」、子の頂点として「根を除いたヒープの要素数(heap_num 列)番目の要素」取り出します。親のデータは、子の数に対応するように複製します。

 推移のアニメーションを作成します。

# ヒープソートのアニメーションを作図
max_h <- floor(log2(N))
d <- 0.1
graph <- ggplot() + 
  geom_path(data = edge_df, 
            mapping = aes(x = coord_x, y = depth, group = edge_id), 
            linewidth = 1) + # 辺
  geom_path(data = heap_edge_df,
            mapping = aes(x = coord_x, y = depth, group = edge_label),
            color = "blue", linewidth = 1, linetype = "dashed", na.rm = TRUE) + # ヒープ対象の辺
  geom_path(data = swap_edge_df,
            mapping = aes(x = coord_x, y = depth, group = operation),
            color = "red", linewidth = 1, linetype = "dashed", na.rm = TRUE) + # 入替対象の辺
  geom_point(data = target_vertex_df,
             mapping = aes(x = coord_x, y = depth, color = vertex_type, group = vertex_label),
             size = 12, alpha = 0.6, na.rm = TRUE) + # 頂点ごとの属性
  geom_point(data = trace_vertex_df,
             mapping = aes(x = coord_x, y = depth, group = id),
             size = 9, shape = "circle filled", fill = "white", stroke = 1) + # 頂点
  geom_text(data = trace_vertex_df,
            mapping = aes(x = coord_x, y = depth, label = as.character(value), group = id),
            size = 3) + # 値ラベル
  geom_text(data = dup_label_df,
            mapping = aes(x = coord_x, y = depth, label = dup_label, group = id),
            size = 3, vjust = 2.6) + # 重複ラベル
  geom_text(data = index_df, 
            mapping = aes(x = coord_x, y = depth, label = index, hjust = offset), 
            size = 3, vjust = -2, color = "green4") + # インデックスラベル
  gganimate::transition_states(states = operation, transition_length = 9, state_length = 1, wrap = FALSE) + # フレーム遷移
  gganimate::ease_aes("cubic-in-out") + # 遷移の緩急
  scale_x_continuous(labels = NULL) + 
  scale_y_reverse(breaks = 0:max_h, minor_breaks = FALSE) + 
  scale_color_manual(breaks = c("max", "heap", "sort"),
                     values = c("red", "blue", "orange"),
                     labels = c("swap", "heap", "sorted"),
                     name = "vertex") + # (凡例表示用)
  coord_cartesian(xlim = c(0, 1), ylim = c(max_h+d, -d)) + 
  labs(title = "heap sort", 
       subtitle = "operation: {next_state}", 
       x = "", y = "depth")

# フレーム数を取得
frame_num <- max(trace_df[["operation"]]) + 1

# 遷移フレーム数を指定
s <- 20

# gif画像を作成
gganimate::animate(
  plot = graph, 
  nframes = (frame_num + 2)*s, start_pause = s, end_pause = s, fps = 20, 
  width = 800, height = 400, 
  renderer = gganimate::gifski_renderer()
)

 ヒープとして扱う要素を青色、ヒープの先頭(最大値)と末尾の要素を赤色、ソート済みの要素をオレンジ色の塗りつぶしで示します。
 最大値の探索(ヒープ化)と入替(削除)を繰り返すことで、降順に(右下から)ソートされるのを確認できます。また、ヒープ化時の入れ替えが、一連の頂点(子孫間)の入れ替えなのが分かります。

 この記事では、ヒープソートを実装しました。次の記事では、バケットソートを実装します。

12.8 バケットソートの実装と可視化

 バケットソート(bucket sort)を実装します。また、ソートの推移のアニメーションを作成して、アルゴリズムを確認します。

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(patchwork)
library(gganimate)

 この記事では、基本的に パッケージ名::関数名() の記法を使っているので、パッケージを読み込む必要はありません。ただし、作図コードについてはパッケージ名を省略しているので、ggplot2 を読み込む必要があります。
 また、ネイティブパイプ演算子 |> を使っています。magrittr パッケージのパイプ演算子 %>% に置き換えても処理できますが、その場合は magrittr を読み込む必要があります。

 「ソートの可視化」において、patchwork パッケージを利用して複数のグラフを並べて描画、gganimate パッケージを利用してアニメーション(gif画像)を作成します。ソートの実装自体には使いません。

ソートの実装

 まずは、バケットソートのアルゴリズムを関数として実装します。

 バケットソートを実装します。

# バケットソートの実装
bucket_sort <- function(vec, max_val = NA) {
  
  # 数列の要素数を取得
  N <- length(vec)
  
  # 最大値を設定
  if(is.na(max_val)) {
    max_val <- max(vec)
  }
  
  # 集計用の中間変数を初期化
  num_vec <- rep(x = 0, times = max_val)
  sum_vec <- rep(x = 0, times = max_val)
  val_vec <- rep(x = NA, times = N)
  
  # 要素ごとに処理
  for(i in 1:N) { # (前から)
    
    # 重複数をカウント
    num_vec[vec[i]] <- num_vec[vec[i]] + 1
  }
  
  # 値ごとに処理
  for(v in 1:max_val) {
    
    # 累積要素数をカウント
    if(v == 1) { # (初回)
      sum_vec[v] <- num_vec[v]
    } else {
      sum_vec[v] <- sum_vec[v-1] + num_vec[v]
    }
  }
  
  # 要素ごとに処理
  for(i in N:1) { # (後から)
    
    # 累積要素数番目に値を格納
    val_vec[sum_vec[vec[i]]] <- vec[i]
    
    # 累積要素数をカウントダウン
    sum_vec[vec[i]] <- sum_vec[vec[i]] - 1
  }
  
  # 数列を更新
  vec <- val_vec
  
  # 数列を出力
  return(vec)
}

 全ての要素が正の整数の数列を vec、数列の要素数を N とします。(本だと0以上の整数ですが、Rはインデックスが1からなので0より大きい整数とします。0を含める場合は、処理の前に数列全体に1を足しておきソート後の出力の前に1を引けばいいと思いますが、この例では扱いません。)
 また、中間変数として扱う際の最大値を max_val として指定します。指定しない場合は、max() で実際の最大値を設定しますが、ソート前に最大値を得るのはおかしいのかもしれません。)

 vec の各要素が正の整数であるのを利用して、vec に格納されている値(1 から max_val までの整数)ごとの要素数をカウントします。初期値が全て 0 のベクトル num_vecvec[i] 番目(i1 から N までの整数)に 1 を足していくことで、要素数が得られます。
 num_vec を用いて累積要素数をカウントします。整数 v-1 までの累積要素数 sum_vec[v-1] と整数 v の要素数 num_vec[v] を足すことで、vec に含まれる要素数の累積和(整数 v までの要素数)が得られます。
 sum_vec を用いて昇順に要素を並べ換えます。i 番目の要素の値(整数) vec[i] の累積要素数 sum_vec[vec[i]] 番目に値 vec[i] を格納し、sum_vec[vec[i]] のカウントを1つ減らしていく(iN から 1 までの(逆順の)整数)ことで、昇順に並んだベクトルを得られます。

 アルゴリズムの詳しい解説は、本のcode 12.5などを参照してください。

 実装した関数を試してみます。

 要素数と最大値を指定して、数列を作成します。

# 要素数を指定
N <- 50

# 最大値を指定
max_val <- 100

# 数列を生成
a <- sample(x = 1:max_val, size = N, replace = TRUE)
a
##  [1] 99 48  2  4 28 20 39 52 33 24 51  2 75  2 22 94 11  1 79 60 95 66 11 60 63
## [26] 38 13 16 58 35 97 11 17  5  7 49 92  5 42 74 90 74 61 22  4 82 77 17 70 78

 乱数生成関数 sample() などで数値ベクトルを生成します。

 バケットソートにより並べ替えます。

# ソート
bucket_sort(a, max_val = max_val)
##  [1]  1  2  2  2  4  4  5  5  7 11 11 11 13 16 17 17 20 22 22 24 28 33 35 38 39
## [26] 42 48 49 51 52 58 60 60 61 63 66 70 74 74 75 77 78 79 82 90 92 94 95 97 99

 目で確認するのはアレなので、組み込み関数のソート結果と比較して確認します。

# 確認
sum(!(bucket_sort(a, max_val) == sort(a)))
## [1] 0

 !TRUE, FALSE を反転させて sum() で合計することで、一致しない要素数を得られます。結果が 0 であれば全ての要素が一致したことが分かります。

ソートの可視化

 次は、バケットソートのアルゴリズムをグラフで確認します。

 数列の初期値を作成して、グラフで確認します。

 要素数と最大値を指定して、数列を作成します。

# 要素数を指定
N <- 30

# 最大値を指定
max_val <- 20

# 数列を生成
a <- sample(x = 1:max_val, size = N, replace = TRUE)
table(a)
## a
##  1  3  4  5  6  7  8  9 11 13 14 16 17 18 19 20 
##  3  2  1  2  1  2  1  2  2  2  1  3  3  1  2  2

 「ソートの実装」のときのコードで数列を作成します。

 数列をデータフレームに格納します。

# 要素数を取得
N <- length(a)

# 数列を格納
seq_df <- tibble::tibble(
  index = 1:N, # 各試行のインデックス
  value = a    # 要素
  #value = bucket_sort(a) # 要素
)
seq_df
## # A tibble: 30 × 2
##    index value
##    <int> <int>
##  1     1     1
##  2     2    19
##  3     3     9
##  4     4    14
##  5     5    17
##  6     6    18
##  7     7    13
##  8     8    17
##  9     9     1
## 10    10    19
## # ℹ 20 more rows

 作図用に、インデックスとあわせて数列を格納します。数列のみが与えられている場合は、要素数を N とします。

 数列を棒グラフで確認します。

# 数列を作図
ggplot() + 
  geom_bar(data = seq_df, 
           mapping = aes(x = index, y = value, fill = factor(value)), stat = "identity") + 
  theme(panel.grid.minor.x = element_blank()) + # 図の体裁
  labs(title = "numerical sequence", 
       fill = "value", 
       x = "index", y = "value")

 この要素(バー)を昇順に並べ替えます。

 簡易版と詳細版の2つのアニメーションでアルゴリズムを確認します。

操作の確認

 作図用に、アルゴリズムに従う途中経過のデータを作成します。

 数列の初期値の描画用のデータフレームを作成します。

# 数列を格納
sequence_df <- tibble::tibble(
  iteration = 0, # フレーム番号
  id        = 1:N, # 要素ID
  # 数列用
  index = 1:N, 
  value = a, 
  # タイル用
  x        = index, 
  y        = 0.5 * value, 
  height   = value, 
  alpha    = 1, 
  linetype = "blank", 
  # ラベル用
  label_y = 0
) |> 
  dplyr::group_by(value) |> # IDの割当用
  dplyr::mutate(
    # 重複ラベル用
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup()
sequence_df |> 
  dplyr::select(!c(height, alpha, linetype, dup_id, dup_num)) # 資料作成用に間引き
## # A tibble: 30 × 8
##    iteration    id index value     x     y label_y dup_label
##        <dbl> <int> <int> <int> <int> <dbl>   <dbl> <chr>    
##  1         0     1     1     1     1   0.5       0 "(1)"    
##  2         0     2     2    19     2   9.5       0 "(1)"    
##  3         0     3     3     9     3   4.5       0 "(1)"    
##  4         0     4     4    14     4   7         0 ""       
##  5         0     5     5    17     5   8.5       0 "(1)"    
##  6         0     6     6    18     6   9         0 ""       
##  7         0     7     7    13     7   6.5       0 "(1)"    
##  8         0     8     8    17     8   8.5       0 "(2)"    
##  9         0     9     9     1     9   0.5       0 "(2)"    
## 10         0    10    10    19    10   9.5       0 "(2)"    
## # ℹ 20 more rows

 数列のインデックスと値を格納して、作図用の値(列)を作成します。
 重複している値の要素にのみ、それぞれ通し番号のラベルを作成します。
 アニメーション用のフレーム番号を 0 とし、以降は通し番号を指定します。

 数列の初期値の棒グラフを作成します。

# 数列を作図
sequence_graph <- ggplot() + 
  geom_tile(data = sequence_df,
            mapping = aes(x = x, y = y, width = 1, height = height, fill = factor(value))) + # 全要素
  # geom_bar(data = sequence_df, 
  #          mapping = aes(x = x, y = value, fill = factor(value)), stat = "identity") + # 全要素
  geom_text(data = sequence_df, 
            mapping = aes(x = x, y = label_y, label = as.character(value)), 
            vjust = -1.5, size = 5) + # 要素ラベル
  geom_text(data = sequence_df, 
            mapping = aes(x = x, y = label_y, label = dup_label), 
            vjust = -0.5, size = 4) + # 重複ラベル
  coord_equal(ratio = 1, xlim = c(1, N), ylim = c(0, max_val)) + # アスペクト比
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "numerical sequence", 
       subtitle = "initial value", 
       fill = "value", 
       x = "index", y = "value")
sequence_graph

 アニメーションの作図の確認として geom_tile() でバーを描画しています。x軸(x 引数)を数列のインデックス、y軸(y 引数)を数列の半分の値、高さ(height 引数)を数列の値にすると、数列に対応した棒グラフを描画できます。
 geom_bar() でバーを描画する場合は、y 引数に value 列を指定します。

 値ごとの要素数の描画用のデータフレームを作成します。

# 値ごとの要素数をカウント
count_df <- tibble::tibble(
  iteration = 1, # フレーム番号
  id        = 1:N, # 要素ID
  # 数列用
  value = a, 
  # タイル用
  x       = value, 
  height  = 1, 
  alpha   = 0.1, 
  linetype = "dashed"
) |> 
  dplyr::group_by(value) |> # IDの割当・カウント用
  dplyr::mutate(
    count = dplyr::row_number(id), # 値ごとの累積要素数
    y     = count - 0.5, 
    # ラベル用
    label_y = y - 0.5, 
    # 重複ラベル用
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup()
count_df |> 
  dplyr::select(!c(height, alpha, linetype, dup_id, dup_num)) # 資料作成用に間引き
## # A tibble: 30 × 8
##    iteration    id value     x count     y label_y dup_label
##        <dbl> <int> <int> <int> <int> <dbl>   <dbl> <chr>    
##  1         1     1     1     1     1   0.5       0 "(1)"    
##  2         1     2    19    19     1   0.5       0 "(1)"    
##  3         1     3     9     9     1   0.5       0 "(1)"    
##  4         1     4    14    14     1   0.5       0 ""       
##  5         1     5    17    17     1   0.5       0 "(1)"    
##  6         1     6    18    18     1   0.5       0 ""       
##  7         1     7    13    13     1   0.5       0 "(1)"    
##  8         1     8    17    17     2   1.5       1 "(2)"    
##  9         1     9     1     1     2   1.5       1 "(2)"    
## 10         1    10    19    19     2   1.5       1 "(2)"    
## # ℹ 20 more rows

 値ごとに(value 列でグループ化して)、row_number() を使って要素IDの小さい順に累積要素数を割り当てます(重複IDの割り当てと同じ処理です)。

 数列に含まれる値ごとの要素数の積み上げ棒グラフを作成します。

# 最大値を取得
max_cnt <- max(count_df[["count"]]) + 1

# 値ごとの要素数を作図
count_graph <- ggplot() + 
  geom_tile(data = count_df,
            mapping = aes(x = x, y = y, width = 1, height = height,
                          fill = factor(value), color = factor(value)),
            alpha = 0.1, linewidth = 0.6, linetype = "dashed") + # 全要素
  # geom_bar(data = count_df, 
  #          mapping = aes(x = x, y = 1, fill = factor(value), color = factor(value)), 
  #          stat = "identity", position = "stack", width = 1, 
  #          alpha = 0.1, linewidth = 0.6, linetype = "dashed") + # 全要素
  geom_text(data = count_df, 
            mapping = aes(x = x, y = label_y, label = as.character(value)), 
            vjust = -1.5, size = 3) + # 要素ラベル
  geom_text(data = count_df, 
            mapping = aes(x = x, y = label_y, label = dup_label), 
            vjust = -0.5, size = 2) + # 重複ラベル
  coord_equal(ratio = 1, xlim = c(1, N), ylim = c(0, max_cnt)) + # アスペクト比
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "count", 
       fill = "value", 
       x = "value", y = "count")
count_graph

 x軸(x 引数)を数列の値、y軸(y 引数)を値ごとの累積要素数から半タイル分ズラした(0.5 を引いた)値、高さ(height 引数)を 1 にすると、数列の要素数に対応した棒グラフを描画できます。
 geom_bar() を使う場合は、y 引数に 1position = "stack" (デフォルト)を指定して積み上げ棒グラフとして描画します。

 累積要素数とインデックスの対応関係の描画用のデータフレームを作成します。

# 累積要素数をカウント
reorder_df <- tibble::tibble(
  iteration = 2, # フレーム番号
  id        = 1:N, # 要素ID
  # 数列用
  value = a, 
  # タイル用
  y        = 0.5, 
  height   = 1, 
  alpha    = 0.1, 
  linetype = "dashed", 
  # ラベル用
  label_y = y - 0.5
) |> 
  dplyr::arrange(value, id) |> # 累積カウント用
  dplyr::mutate(
    x = dplyr::row_number() # 累積要素数
  ) |> 
  dplyr::arrange(id) |> 
  dplyr::group_by(value) |> # IDの割当用
  dplyr::mutate(
    # 重複ラベル用
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup()
reorder_df |> 
  dplyr::select(!c(height, alpha, linetype, dup_id, dup_num)) # 資料作成用に間引き
## # A tibble: 30 × 7
##    iteration    id value     y label_y     x dup_label
##        <dbl> <int> <int> <dbl>   <dbl> <int> <chr>    
##  1         2     1     1   0.5       0     1 "(1)"    
##  2         2     2    19   0.5       0    27 "(1)"    
##  3         2     3     9   0.5       0    13 "(1)"    
##  4         2     4    14   0.5       0    19 ""       
##  5         2     5    17   0.5       0    23 "(1)"    
##  6         2     6    18   0.5       0    26 ""       
##  7         2     7    13   0.5       0    17 "(1)"    
##  8         2     8    17   0.5       0    24 "(2)"    
##  9         2     9     1   0.5       0     2 "(2)"    
## 10         2    10    19   0.5       0    28 "(2)"    
## # ℹ 20 more rows

 値と要素IDの小さい順に、row_number() を使って通し番号を割り当てます。

 昇順の値ごとの累積要素数の棒グラフを作成します。

# 累積要素数を作図
reorder_graph <- ggplot() + 
  geom_tile(data = reorder_df,
            mapping = aes(x = x, y = y, width = 1, height = height,
                          fill = factor(value), color = factor(value)),
            alpha = 0.1, linewidth = 0.6, linetype = "dashed") + # 全要素
  # geom_bar(data = reorder_df, 
  #          mapping = aes(x = x, y = 1, fill = factor(value), color = factor(value)), 
  #          stat = "identity", width = 1, 
  #          alpha = 0.1, linewidth = 0.6, linetype = "dashed") + # 全要素
  geom_text(data = reorder_df, 
            mapping = aes(x = x, y = label_y, label = as.character(value)), 
            vjust = -1.5, size = 3) + # 要素ラベル
  geom_text(data = reorder_df, 
            mapping = aes(x = x, y = label_y, label = dup_label), 
            vjust = -0.5, size = 2) + # 重複ラベル
  coord_equal(ratio = 1, xlim = c(1, N), ylim = c(0, max_cnt)) + # アスペクト比
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "cumulative sum", 
       fill = "value", 
       x = "value", y = "count")
reorder_graph

 x軸(x 引数)を累積要素数(昇順の通し番号)、y軸(y 引数)を 0.5、高さ(height 引数)を 1 にすると、昇順のインデックスに対応した棒グラフを描画できます。
 geom_bar() を使う場合は、y 引数に 1 を指定します。

 昇順にソートした数列の描画用のデータフレームを作成します。

# 数列をソート
swap_df <- tibble::tibble(
  iteration = 3, # フレーム番号
  id        = 1:N, # 要素ID
  # 数列用
  value = a, 
  # タイル用
  y        = 0.5 * value, 
  height   = value, 
  alpha    = 1, 
  linetype = "blank", 
  # ラベル用
  label_y = 0
) |> 
  dplyr::arrange(value, id) |> # 入替用
  dplyr::mutate(
    index = dplyr::row_number(), # 入替後のインデックス
    x     = index
  ) |> 
  dplyr::arrange(id) |> 
  dplyr::group_by(value) |> # IDの割当用
  dplyr::mutate(
    # 重複ラベル用
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup()
swap_df |> 
  dplyr::select(!c(height, alpha, linetype, dup_id, dup_num)) # 資料作成用に間引き
## # A tibble: 30 × 8
##    iteration    id value     y label_y index     x dup_label
##        <dbl> <int> <int> <dbl>   <dbl> <int> <int> <chr>    
##  1         3     1     1   0.5       0     1     1 "(1)"    
##  2         3     2    19   9.5       0    27    27 "(1)"    
##  3         3     3     9   4.5       0    13    13 "(1)"    
##  4         3     4    14   7         0    19    19 ""       
##  5         3     5    17   8.5       0    23    23 "(1)"    
##  6         3     6    18   9         0    26    26 ""       
##  7         3     7    13   6.5       0    17    17 "(1)"    
##  8         3     8    17   8.5       0    24    24 "(2)"    
##  9         3     9     1   0.5       0     2     2 "(2)"    
## 10         3    10    19   9.5       0    28    28 "(2)"    
## # ℹ 20 more rows

 値と要素IDの小さい順に、row_number() を使ってインデックス(通し番号)を割り当てます。

 昇順に並べ替えた数列の棒グラフを作成します。

# ソート済み数列を作図
swap_graph <- ggplot() + 
  geom_tile(data = swap_df,
            mapping = aes(x = x, y = y, width = 1, height = height, fill = factor(value))) + # 全要素
  # geom_bar(data = swap_df, 
  #          mapping = aes(x = x, y = value, fill = factor(value)), 
  #          stat = "identity") + # 全要素
  geom_text(data = swap_df, 
            mapping = aes(x = x, y = label_y, label = as.character(value)), 
            vjust = -1.5, size = 5) + # 要素ラベル
  geom_text(data = swap_df, 
            mapping = aes(x = x, y = label_y, label = dup_label), 
            vjust = -0.5, size = 4) + # 重複ラベル
  coord_equal(ratio = 1, xlim = c(1, N), ylim = c(0, max_val)) + # アスペクト比
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "bucket sort", 
       fill = "value", 
       x = "index", y = "value")
swap_graph

 x軸(x 引数)を累積要素数(昇順の通し番号)、y軸(y 引数)を数列の半分の値、高さ(height 引数)を数列の値にすると、数列に対応した棒グラフを描画できます。
 geom_bar() でバーを描画する場合は、y 引数に value 列を指定します。

 各操作のグラフを並べて推移を確認します。

 4つのグラフを並べて描画します。

# 並べて描画
patchwork::wrap_plots(
  sequence_graph, swap_graph, 
  count_graph, reorder_graph, 
  nrow = 2, ncol = 2, widths = 1, heights = c(1, max_cnt/max_val)
)

 wrap_plots() を使って複数のグラフを1つにまとめて描画します。

 各操作のグラフを切り替えて推移を確認します。

 4つのデータフレームをまとめます。

# データを結合
trace_df <- dplyr::bind_rows(
  sequence_df, 
  count_df, 
  reorder_df, 
  swap_df
)
trace_df |> 
  dplyr::arrange(iteration, id)
## # A tibble: 120 × 14
##    iteration    id index value     x     y height alpha linetype label_y dup_id
##        <dbl> <int> <int> <int> <int> <dbl>  <dbl> <dbl> <chr>      <dbl>  <int>
##  1         0     1     1     1     1   0.5      1     1 blank          0      1
##  2         0     2     2    19     2   9.5     19     1 blank          0      1
##  3         0     3     3     9     3   4.5      9     1 blank          0      1
##  4         0     4     4    14     4   7       14     1 blank          0      1
##  5         0     5     5    17     5   8.5     17     1 blank          0      1
##  6         0     6     6    18     6   9       18     1 blank          0      1
##  7         0     7     7    13     7   6.5     13     1 blank          0      1
##  8         0     8     8    17     8   8.5     17     1 blank          0      2
##  9         0     9     9     1     9   0.5      1     1 blank          0      2
## 10         0    10    10    19    10   9.5     19     1 blank          0      2
## # ℹ 110 more rows
## # ℹ 3 more variables: dup_num <int>, dup_label <chr>, count <int>

 bind_rows() で複数のデータフレームを行方向に結合します。

 4つの操作のアニメーションを作成します。

# ソートのアニメーションを作図
graph <- ggplot() + 
  geom_tile(data = trace_df, 
            mapping = aes(x = x, y = y, width = 1, height = height, 
                          fill = factor(value), color = factor(value), alpha = alpha, 
                          linetype = linetype, group = factor(id)), 
            linewidth = 1) + # 全要素
  geom_text(data = trace_df, 
            mapping = aes(x = x, y = label_y, label = as.character(value), group = factor(id)), 
            vjust = -1.5, size = 5) + # 要素ラベル
  geom_text(data = trace_df, 
            mapping = aes(x = x, y = label_y, label = dup_label, group = factor(id)), 
            vjust = -0.5, size = 4) + # 重複ラベル
  gganimate::transition_states(states = iteration, transition_length = 9, state_length = 1, wrap = FALSE) + # フレーム遷移
  gganimate::ease_aes("cubic-in-out") + # 遷移の緩急
  scale_linetype_identity() + 
  scale_linewidth_identity() + 
  scale_y_continuous(sec.axis = sec_axis(trans = ~-., name = "count")) + # (操作の確認用)
  coord_equal(ratio = 1) + # アスペクト比
  theme(panel.grid.minor.x = element_blank(), 
        legend.position = "none") + 
  labs(title = "bucket sort", 
       subtitle = "iteration: {next_state}", 
       fill = "value", 
       x = "index, value", y = "value, count")

# フレーム数を取得
frame_num <- trace_df[["iteration"]] |> 
  (\(vec) {max(vec) + 1})()

# 遷移フレーム数を指定
s <- 20

# gif画像を作成
gganimate::animate(
  plot = graph, 
  nframes = (frame_num + 2)*s, start_pause = s, end_pause = s, fps = 20, 
  width = 1000, height = 800, 
  renderer = gganimate::gifski_renderer()
)

 transition_states() のフレーム制御の引数 states にフレーム番号列 iteration を指定して、要素の値や要素数の対応関係が補完されるように geom_tile() で棒グラフを描画します。
 animate()plot 引数にグラフオブジェクト、nframes 引数にフレーム数を指定して、gif画像を作成します。
 フレーム数は 4 です。transition_states() でアニメーションを作成すると、フレーム間の値を補完して、状態の遷移を描画します。遷移のフレームを s として、試行回数の s 倍の値を nframes 引数に指定します。

 元の「数列の値」、「値ごとの要素数」、昇順にカウントした「累積要素数」、ソート後の「数列のインデックス」の対応関係を確認できます。

入替の確認

 続いて「ソートの実装」のアルゴリズムに対応するように、数列の前の要素からカウントし、後の要素から配置していくデータを作成します。

 カウント時の数列の描画用のデータフレームを作成します。

# インデックスの小さい順に要素を取り出す
sequence_df <- tidyr::expand_grid(
  iteration = 0:N, # カウント時の試行回数
  id        = 1:N  # 元のインデックス
) |> # 試行ごとに要素IDを複製
  dplyr::mutate(
    # 数列用
    index = id, 
    value = a[id], 
    # タイル用
    x        = index, 
    y        = 0.5 * value, 
    height   = value, 
    alpha    = 1, 
    linetype = "blank", 
    # ラベル用
    label_y = 0
  ) |> 
  dplyr::group_by(iteration, value) |> # IDの割当用
  dplyr::mutate(
    # 重複ラベル用
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup() |> 
  dplyr::filter(iteration < id) # 未カウントの要素を抽出
sequence_df |> 
  dplyr::select(!c(height, alpha, linetype, dup_id, dup_num)) # 資料作成用に間引き
## # A tibble: 465 × 8
##    iteration    id index value     x     y label_y dup_label
##        <int> <int> <int> <int> <int> <dbl>   <dbl> <chr>    
##  1         0     1     1     1     1   0.5       0 "(1)"    
##  2         0     2     2    19     2   9.5       0 "(1)"    
##  3         0     3     3     9     3   4.5       0 "(1)"    
##  4         0     4     4    14     4   7         0 ""       
##  5         0     5     5    17     5   8.5       0 "(1)"    
##  6         0     6     6    18     6   9         0 ""       
##  7         0     7     7    13     7   6.5       0 "(1)"    
##  8         0     8     8    17     8   8.5       0 "(2)"    
##  9         0     9     9     1     9   0.5       0 "(2)"    
## 10         0    10    10    19    10   9.5       0 "(2)"    
## # ℹ 455 more rows

 試行番号(初期値を 0 として N までの整数)と要素ID(初期値のインデックス)(1 から N までの整数)の全ての組み合わせを expand_grid() で作成します。
 「操作の確認」のときと同様に処理して、未カウントの要素(試行番号よりIDが大きい行)を抽出します。

 カウント時の値ごとの要素数の描画用のデータフレームを作成します。

# インデックスの小さい順に要素を取り出してカウント
count_df <- tidyr::expand_grid(
  iteration = 0:N, # カウント時の試行回数
  id        = 1:N  # 元のインデックス
) |> # 試行ごとに要素IDを複製
  dplyr::mutate(
    # 数列用
    value = a[id], 
    # タイル用
    x        = value, 
    height   = 1, 
    alpha    = 0.1, 
    linetype = "dashed"
  ) |> 
  dplyr::group_by(iteration, value) |> # IDの割当・カウント用
  dplyr::mutate(
    # タイル用
    count = dplyr::row_number(id), # 値ごとの要素数
    y     = -count + 0.5, 
    # ラベル用
    label_y = y - 0.5, 
    # 重複ラベル用
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup() |> 
  dplyr::filter(iteration >= id) # カウント済みの要素を抽出
count_df |> 
  dplyr::select(!c(height, alpha, linetype, dup_id, dup_num)) # 資料作成用に間引き
## # A tibble: 465 × 8
##    iteration    id value     x count     y label_y dup_label
##        <int> <int> <int> <int> <int> <dbl>   <dbl> <chr>    
##  1         1     1     1     1     1  -0.5      -1 "(1)"    
##  2         2     1     1     1     1  -0.5      -1 "(1)"    
##  3         2     2    19    19     1  -0.5      -1 "(1)"    
##  4         3     1     1     1     1  -0.5      -1 "(1)"    
##  5         3     2    19    19     1  -0.5      -1 "(1)"    
##  6         3     3     9     9     1  -0.5      -1 "(1)"    
##  7         4     1     1     1     1  -0.5      -1 "(1)"    
##  8         4     2    19    19     1  -0.5      -1 "(1)"    
##  9         4     3     9     9     1  -0.5      -1 "(1)"    
## 10         4     4    14    14     1  -0.5      -1 ""       
## # ℹ 455 more rows

 「操作の確認」のときと同様に処理して、カウント済みの要素(試行番号がID以上の行)を抽出します。ただし、バーが重ならないようにカウントの符号を反転させた(-1 を掛けた)値をy軸の値とします。

 昇順に配置時の累積要素数とインデックスの対応関係の描画用のデータフレームを作成します。

# インデックスの小さい順に要素を戻してカウント
discount_df <- tidyr::expand_grid(
  iteration = 0:N + N+1, # 入替時の試行回数
  id        = 1:N # 元のインデックス
) |> # 試行ごとに要素IDを複製
  dplyr::mutate(
    # 数列用
    value = a[id], 
    # タイル用
    y        = -0.5, 
    height   = 1, 
    alpha    = 0.1, 
    linetype = "dashed", 
    # ラベル用
    label_y = y - 0.5
  ) |> 
  dplyr::arrange(iteration, value, id) |> # 累積カウント用
  dplyr::group_by(iteration) |> # 累積カウント用
  dplyr::mutate(
    # タイル用
    x = dplyr::row_number() # 累積要素数
  ) |> 
  dplyr::arrange(iteration, id) |> 
  dplyr::group_by(iteration, value) |> # IDの割当用
  dplyr::mutate(
    # 重複ラベル用
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup() |> 
  dplyr::filter(iteration-(N+1) < N+1-id) # 未入替の要素を抽出
discount_df |> 
  dplyr::select(!c(height, alpha, linetype, dup_id, dup_num)) # 資料作成用に間引き
## # A tibble: 465 × 7
##    iteration    id value     y label_y     x dup_label
##        <dbl> <int> <int> <dbl>   <dbl> <int> <chr>    
##  1        31     1     1  -0.5      -1     1 "(1)"    
##  2        31     2    19  -0.5      -1    27 "(1)"    
##  3        31     3     9  -0.5      -1    13 "(1)"    
##  4        31     4    14  -0.5      -1    19 ""       
##  5        31     5    17  -0.5      -1    23 "(1)"    
##  6        31     6    18  -0.5      -1    26 ""       
##  7        31     7    13  -0.5      -1    17 "(1)"    
##  8        31     8    17  -0.5      -1    24 "(2)"    
##  9        31     9     1  -0.5      -1     2 "(2)"    
## 10        31    10    19  -0.5      -1    28 "(2)"    
## # ℹ 455 more rows

 カウント時に続く試行番号(N+1 から 2*N+1 までの整数)と要素ID(1 から N までの整数)の全ての組み合わせを expand_grid() で作成します。
 「操作の確認」のときと同様に処理して、入替前の要素(「試行番号」より「後から数えたID」が大きい行)を抽出します。

 昇順に配置時の数列の描画用のデータフレームを作成します。

# インデックスの小さい順に要素を入れ替えて戻す
swap_df <- tidyr::expand_grid(
  iteration = 0:N + N+1, # 入替時の試行回数
  id        = 1:N # 元のインデックス
) |> # 試行ごとに要素IDを複製
  dplyr::mutate(
    # 数列用
    value = a[id], 
    # タイル用
    y        = 0.5 * value, 
    height   = value, 
    alpha    = 1, 
    linetype = "blank", 
    # ラベル用
    label_y = 0
  ) |> 
  dplyr::arrange(iteration, value, id) |> # 入替用
  dplyr::group_by(iteration) |> # 入替用
  dplyr::mutate(
    # タイル用
    index = dplyr::row_number(), # 入替後のインデックス
    x     = index
  ) |> 
  dplyr::arrange(iteration, id) |> 
  dplyr::group_by(iteration, value) |> # IDの割当用
  dplyr::mutate(
    # 重複ラベル用
    dup_id    = dplyr::row_number(id), # 重複IDを割り当て
    dup_num   = max(dup_id), # 重複の判定用
    dup_label = dplyr::if_else(
      condition = dup_num > 1, true = paste0("(", dup_id, ")"), false = ""
    ) # 重複要素のみラベルを作成
  ) |> 
  dplyr::ungroup() |> 
  dplyr::filter(iteration-(N+1) >= N+1-id) # 入替済みの要素を抽出
swap_df |> 
  dplyr::select(!c(height, alpha, linetype, dup_id, dup_num)) # 資料作成用に間引き
## # A tibble: 465 × 8
##    iteration    id value     y label_y index     x dup_label
##        <dbl> <int> <int> <dbl>   <dbl> <int> <int> <chr>    
##  1        32    30     5   2.5       0     8     8 "(2)"    
##  2        33    29     8   4         0    12    12 ""       
##  3        33    30     5   2.5       0     8     8 "(2)"    
##  4        34    28    20  10         0    30    30 "(2)"    
##  5        34    29     8   4         0    12    12 ""       
##  6        34    30     5   2.5       0     8     8 "(2)"    
##  7        35    27    11   5.5       0    16    16 "(2)"    
##  8        35    28    20  10         0    30    30 "(2)"    
##  9        35    29     8   4         0    12    12 ""       
## 10        35    30     5   2.5       0     8     8 "(2)"    
## # ℹ 455 more rows

 「操作の確認」のときと同様に処理して、入替前の要素(「試行番号」が「後から数えたID」以上の行)を抽出します。

 4つのデータフレームをまとめます。

# データを結合
trace_df <- dplyr::bind_rows(
  sequence_df, 
  count_df, 
  discount_df, 
  swap_df
)
trace_df
## # A tibble: 1,860 × 14
##    iteration    id index value     x     y height alpha linetype label_y dup_id
##        <dbl> <int> <int> <int> <int> <dbl>  <dbl> <dbl> <chr>      <dbl>  <int>
##  1         0     1     1     1     1   0.5      1     1 blank          0      1
##  2         0     2     2    19     2   9.5     19     1 blank          0      1
##  3         0     3     3     9     3   4.5      9     1 blank          0      1
##  4         0     4     4    14     4   7       14     1 blank          0      1
##  5         0     5     5    17     5   8.5     17     1 blank          0      1
##  6         0     6     6    18     6   9       18     1 blank          0      1
##  7         0     7     7    13     7   6.5     13     1 blank          0      1
##  8         0     8     8    17     8   8.5     17     1 blank          0      2
##  9         0     9     9     1     9   0.5      1     1 blank          0      2
## 10         0    10    10    19    10   9.5     19     1 blank          0      2
## # ℹ 1,850 more rows
## # ℹ 3 more variables: dup_num <int>, dup_label <chr>, count <int>


 カウント時と配置時のグラフを切り替えて推移を確認します。

 「操作の確認」のときのコードで作図します。「初期値」の表示と「値ごとの要素数と累積要素数によるインデックス」の遷移を含めるため、フレーム数は 2*N + 2 です。
 バーが重ならないように、要素数のカウントを逆向きの積み上げ棒グラフで表しています。

 左(1番目)の要素から順番にカウントしていき、右(N番目)の要素から逆順に累積要素数に応じて配置していくことでソートされるのを確認できます。また、重複要素の順序が入れ替わらない(安定性・stable性)も確認できます。

 ここでの試行回数(iterationの値)は、この可視化方法における試行番号(作図処理上の値)であり他の手法と対応していません。

 この記事では、バケットソートを実装しました。

参考文献

  • 大槻兼資(著), 秋葉拓哉(監修)『問題解決力を鍛える! アルゴリズムとデータ構造』講談社サイエンティク, 2021年.